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Preface 


The  numerical  simulations  in  this  report  were  conducted  using  Department  of 
Defense  (DoD)  supercomputer  time  provided  through  several  Challenge  Allocations 
awarded  by  the  High  Performance  Computing  Modernization  Office  (HPCMO).  Systems 
at  the  U.S.  Engineer  Research  and  Development  Center  (ERDC)  and  the  Naval 
Oceanographic  Office  (NAVO)  were  used.  The  rawinsonde  data  analyzed  was  obtained 
through  our  participation  in  the  CASES-99  and  VTMX  field  programs,  sponsored  by  the 
National  Science  Foundation  (NSF)  and  the  Department  of  Energy  (DOE).  Additional 
Air  Force  rawinsonde  data  were  provided  by  Dr.  Frank  Ruggiero  of  the  Air  Force 
Research  Laboratory  (AFRL).  The  comparisons  with  aircraft  data  were  made  possible  by 
results  and  analysis  provided  by  Dr.  Donald  Wroblewski  of  Boston  University  (BU). 
Some  of  the  simulations  and  visualization  work  was  done  by  Drs.  Ling  Wang,  Kam  Wan, 
Chris  Bizon,  Chris  Meyer,  and  Ms.  Melissa  Davis-Mansour. 


1.  INTRODUCTION 


In  the  summer  of  2002,  Northwest  Research  Associates,  Inc.  (NWRA)  was  awarded 
a  three-year  contract  by  the  Air  Force  Research  Laboratory  (AFRL)  for  a  proposal 
submitted  to  the  High-Energy  Laser  -  Joint  Tactical  Office  (HEL-JTO)  program  to 
perform  the  research  and  development  needed  to  develop  and  improve  reliable 
Atmospheric  Decision  Aids  (AD As)  for  Air  Force  optical  turbulence  forecasting 
applications  and  operational  support.  Example  applications  include  laser-weapons  and 
laser-communication  systems. 

Though  the  HEL-JTO  program  was  unexpectedly  tenninated  less  than  one  year  after 
project  initiation,  AFRL  continued  to  sponsor  this  effort  under  the  original  contract, 
extending  the  period  of  performance  by  2.5  years  and  fully  funding  the  original  contract 
amount.  We  are  grateful  that  AFRL  managed  to  keep  this  effort  alive  which  made  the 
research  presented  here  possible. 

During  this  project,  we  completed  a  series  of  Direct-Numerical  and  Large-Eddy 
Simulations  (DNS  and  LES)  of  relevant  atmospheric  dynamical  processes  occurring  in 
the  troposphere  and  stratosphere  (e.g.,  wind-shear  instability  and  gravity  wave  breaking), 
and  analyzed  the  statistics  of  the  resulting  atmospheric  turbulence.  Two  DoD  High 
Performance  Computing  Modernization  Program  (HPCMP)  Challenge  Allocations  at  two 
Major  Shared  Resource  Centers  (MSRCs),  i.e.,  ERDC  and  NAVO,  were  awarded  to  us  to 
carry  out  this  part  of  the  work.  We  also  analyzed  hundreds  of  rawinsonde  profiles  in 
order  to  determine  the  parameters  associated  with  turbulence  layers  observed  in  the 
troposphere  and  stratosphere.  This  data  was  obtained  through  our  participation  in  the 
CASES-99  and  VTMX  field  campaigns,  as  well  as  Air  Force  rawinsonde  data  provided 
by  Dr.  Frank  Ruggiero.  Finally,  we  developed  a  Bayesian  Hierarchical  Modeling  (BHM) 
framework  for  combining  the  simulation  and  observational  results  so  that  they  may  be 
used  to  develop  and/or  improve  atmospheric  turbulence  and  optical  turbulence  AD  As. 
This  approach  and  data  can  also  be  used  to  develop  subgrid-scale  turbulence  and  optical 
turbulence  parameterizations  for  operational  mesoscale  Numerical  Weather  Prediction 
(NWP)  models,  e.g.,  the  Weather  Research  and  Forecasting  Model  (WRF). 


2.  BACKGROUND 


At  the  outset  of  this  project,  the  Air  Force  was  developing  an  ADA  to  provide  optical 
turbulence  predictions  in  support  of  the  Airborne  Laser  (ABL).  The  starting  point  for  the 
ABL  ADA  was  the  following  model  of  the  index-of-refraction  second-order  structure 
constant  Cn2: 

(i) 

\kmJ 

where  a=2.8  is  a  universal  constant;  Km  and  Kn  are  eddy  diffusion  coefficients  for 
momentum  and  the  index  of  refraction  n;  <3Z  n  is  the  mean  vertical  gradient  of  n;  and  £0  is 
the  “outer  length  scale”  for  turbulent  motion,  above  which  effects  of  stratification  are  felt. 
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Dewan  (1980)  and  others  argue  that  tQ  is  of  order  and  £0  «  0. 1  L,  where  L  is  the  depth  of  a 

turbulent  mixing  layer.  Equation  (1)  was  proposed  by  Tatarski  in  1961  (Tatarski  1961);  a 
detailed  description  of  its  derivation  is  given  in  Dewan  (1980).  The  critical  three 
assumptions  made  when  deriving  (1)  are:  1)  stationary  turbulence,  in  which  production 
balances  dissipation  (i.e.,  P  ~  s)  for  both  n  and  momentum;  2)  homogeneous  turbulence, 
for  which  all  turbulent  diffusion  terms  play  a  minor  role  and  do  not  influence  the 
dominant  balance  of  P  «  s;  and  3)  Prandtl  mixing-length  arguments,  for  which  transport 
and  mixing  are  described  by  an  eddy  viscosity  and  down-gradient  fluxes.  All  of  the 
quantities  in  (1)  can  be  either  estimated  or  deduced  from  mesoscale-model  output,  except 
for  £0,  which  Dewan  et  al.  (1993)  estimate  from  filtered  wind  profiles.  They  use  a  filter 

width  of  300m,  i.e.,  the  same  resolution  as  archived  radiosonde  data  and  also  the 
resolution  of  the  Air  Force  Weather  Agency  (AFWA)  mesoscale  model  at  ABF 
operational  altitudes,  to  produce  a  linear  fit  of  Y  versus  Sraw,  where  Sraw  is  the  filtered 
shear  rate  and  10Y=(F4/3).  The  large  scatter  in  the  data  (see  Figures  2  and  3  of  Dewan  et 
al.  1993)  indicates  that  F  correlates  poorly  with  the  filtered  shear  rate  Sraw,  and  that  the 
uncertainty  in  Y  is  as  large  as  ±1.  This  translates  to  an  uncertainty  in  Cn  that  spans  as 
much  as  a  factor  of  100,  which  is  as  large  as  (and  in  some  cases  larger  than!)  the  natural 
variation  of  Cn2  exhibited  by  thermosonde  data  (Coulman  et  al.  1995).  As  a  result,  (1) 
does  not  predict  Cn  with  any  greater  fidelity  than  a  compiled  Cn  climatology,  e.g., 
CFEAR-1  (Beland  1993). 

The  reasons  for  the  inability  of  (1)  to  provide  greater  predictive  skill  lies  in  a)  the 
coarseness  of  (1)  as  a  model  for  stratified  turbulence  and  b)  the  absence  of  more  detailed 
physics  in  the  regression  relating  F  and  Sraw.  For  example,  local  maxima  in  Cn~  occur  at 
the  edges  of  stratified  mixing  layers  (Gibson-Wilde  et  al.  2000,  Weme  and  Fritts  2001, 
Pettersson-Reif  et  al.  2002).  These  regions  define  the  boundaries  between  vigorous 
turbulence  in  the  middle  of  the  layers  and  the  laminar  motions  outside.  Based  on  recent 
research  (see  §4),  we  know  these  regions  to  be  nonstationary,  inhomogeneous,  and 
anisotropic  (Werne  and  Fritts  2000,  2001).  And  from  simulation  work  done  as  a  part  of 
this  project,  we  also  know  that  turbulent  diffusion,  counter-gradient  fluxes,  and 
anomalous  pressure-strain  correlations  are  all  significant  here  as  well  (Pettersson-Reif  et 
al.  2002,  Weme  et  al.  2005).  Hence,  all  three  of  the  assumptions  used  to  derive  (1)  are 
violated  in  the  regions  that  are  most  important  for  evaluating  Cn  .  It  is  clear  that  a  more 
quantitative  description  of  Cn2  than  (1)  is  required  if  reliable  forecasts  are  to  be  made. 

3.  TECHNICAL  APPROACH 

During  the  course  of  this  work,  we  pursued  three  distinct  project  components:  1) 
high-resolution  numerical  simulations;  2)  rawinsonde  and  aircraft  data  analysis;  and  3) 
the  development  of  a  new  probabilistic  modeling  approach,  which  can  be  used  both  for 
Cn2  modeling  and  as  the  basis  for  a  new  SGS  parameterization  method  for  real-time 
mesoscale  forecasting.  Below,  we  separately  describe  the  technical  approach  used  for 
each  of  these  components. 
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3.1.  High-Resolution  Simulation  Techniques 


The  technical  approach  employed  for  the  high-resolution  simulations  of  atmospheric 
wind  shears  and  gravity  wave  breaking  used  both  DNS  and  LES  techniques.  The  DNS 
code,  named  TRIPLE,  is  an  incompressible  pseudo-spectral  Navier-Stokes  solver  that 
uses  Fourier  and  Chebyshev  spatial  discretizations.  It  is  a  mature  production  code  that 
has  served  as  the  workhorse  for  six  DoD  HPCMP  and  NSF  Grand  Challenge  projects  and 
two  HPCMP  Capability  Applications  Projects  (CAP)  Phase  II  projects  at  seven  different 
supercomputer  centers  (ERDC,  NAVO,  ARL,  ASC,  AHPCRC,  PSC,  and  SDSC)1  on  ten 
different  MPP  platforms  (Cray  T3D,  T3E,  XT3,  IBM  P3,  P4+,  P5+,  SGI  02K,  03K, 
Onyx,  Compaq  SC40/45).  The  numerical  method  advances  the  Fourier  and/or 
Chebyshev  coefficients  in  spectral  space.  Potential  truncation  errors  are  removed  using 
the  “tau”-correction  (Kleisser  and  Schumann  1980,  Werne  1995,  1998),  and  vertical 
stacking  of  Chebyshev  domains  is  possible.  The  method  is  Nth-order  accurate  in  space, 
where  N  is  the  number  of  spectral  modes.  Depending  on  the  application  (i.e.,  wave¬ 
breaking,  shear,  convection,  etc.),  different  spatial  representations  are  used.  Possible 
boundary  conditions  include  stress-free,  no-slip,  fixed-temperature  or  heat-flux,  radiative, 
incident  gravity  waves,  or  any  combination  thereof.  For  shear- instability  and  reduced- 
equation  studies,  greater  computation  speeds  are  achieved  by  employing  simpler 
boundary  conditions  pennitted  by  these  problems.  For  these  cases,  stress-free  top  and 
bottom  boundaries  permit  Fourier  representations  in  the  vertical,  and  complicated 
inversion  operations  can  be  avoided. 

Time-stepping  is  conducted  using  the  mixed  implicit/explicit  3rd-order  Runge-Kutta 
(RK3)  scheme  of  Spalart  et  al.  (1991),  which  has  the  same  storage  requirements  of  most 
2nd-order  schemes,  allowing  us  to  achieve  larger  numerical  domains  than  would 
otherwise  be  possible.  Nonlinear  terms  are  treated  explicitly,  while  buoyancy  and 
diffusion  are  handled  implicitly.  Derivatives  (nonlinear  terms)  are  computed  efficiently 
in  spectral  (physical)  space. 

The  majority  of  the  LES  work  done  on  this  project  was  accomplished  using  an 
updated  version  of  Dr.  Peter  Sullivan’s  Planetary  Boundary  Layer  (PBL)  code.  Dr. 
Sullivan  was  a  consultant  on  this  project.  Sullivan’s  LES  code  has  been  successfully 
applied  to  a  wide  range  of  problems  in  both  atmospheric  and  oceanic  planetary  boundary 
layers.  The  solution  algorithm  is  a  mixed  pseudo-spectral/finite-difference  collocation 
scheme.  Vertical  derivatives  are  approximated  with  2nd-order  centered  differencing. 
Horizontal  velocity  components,  virtual  potential  temperature,  and  pressure  are  stored  at 
co-located  grid  locations  staggered  vertically  from  the  vertical  velocity  and  SGS  energy. 
Time  integration  uses  an  explicit  3rd-order  Runge-Kutta  scheme  for  all  terms  except 
pressure.  A  fractional-step  method  is  used  to  determine  the  pressure  at  each  Runge-Kutta 
substep.  The  pressure  satisfies  a  Poisson  equation,  which  is  solved  using  Fourier  methods 


1  ARL  is  the  Army  Research  Laboratory;  ASC  is  the  Aeronautical  Systems  Center; 
AHPCRC  is  the  Army  High  Perfonnance  Computing  Research  Center;  PSC  is  the 
Pittsburgh  Supercomputing  Center;  SDSC  is  the  San  Diego  Supercomputing  Center. 
ERDC,  NAVO,  and  ASC  are  DoD  Major  Share  Resource  Centers  (MSRCs). 
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and  the  inversion  of  a  tri-diagonal  matrix.  The  code  is  written  in  MPI  and  has  been  in 
production  on  IBM  SP  and  Cray  platforms  at  NAVO  and  ERDC.  For  current  studies, 
SGS  turbulence  is  modeled  using  either  a  two-part  eddy-viscosity  model  (Sullivan  et  al. 
1994)  requiring  the  solution  of  a  transport  equation  for  SGS  energy  or  a  dynamic  SGS 
procedure  which  uses  either  the  two-part  procedure  or  a  traditional  eddy-viscosity  model 
(e.g.,  Smagorinsky  1963)  and  a  gradient-diffusion  eddy-diffusivity  model.  Alternate  SGS 
prescriptions  can  be  accommodated  in  a  straightforward  fashion.  Horizontal  boundary 
conditions  are  periodic.  The  upper  and  lower  boundaries  may  be  periodic,  stress  free,  or 
radiate  gravity  waves.  Other  vertical  boundary  conditions  can  be  implemented  easily.  A 
novel  feature  of  the  present  LES  code  is  the  ability  to  nest  high-resolution  meshes  within 
the  larger  outer  domain  (Sullivan  et  al.  1996).  Two-way  communication  between  nested 
and  outer  grids  utilizes  spectral  interpolation  and  Gennano’s  identity  (Germano  et  al. 
1991),  which  enforces  total-flux  conservation  on  the  different  meshes. 

In  order  to  isolate  the  LES  SGS  model  during  testing  and  comparisons  with  the  DNS 
solutions,  in  the  last  year  of  this  project  we  converted  TRIPLE  from  a  DNS-only  code  to 
a  DNS/LES  code,  which  could  include  the  SGS  schemes  in  Dr.  Sullivan’s  PBL  code  by 
setting  an  LES  parameter  in  TRIPLE’s  job-preparation  input  file. 

3.1.1.  TRIPLE  Optimization  Details 

Three-dimensional  Fast  Fourier  Transforms  (FFTs)  comprise  75%  of  TRIPLE’s  cost, 
and  care  is  therefore  taken  to  ensure  they  are  performed  efficiently.  Operation  counts  are 
reduced  by  exploiting  the  symmetry  properties  of  real  FFTs,  and  loop-level  optimization 
and  temporary  storage  are  used  to  maximize  cache  reuse,  e.g.,  trigonometric  multi-angle 
recurrence  relations  replace  traditional  table-look-up  algorithms,  and  higher  radix 
routines  (e.g.,  powers  of  2,  3,  4,  and  5)  are  separately  coded.  Finally,  transposes  are  used 
to  ensure  FFTs  are  only  conducted  in  the  contiguous  array  direction. 

Inter-processor  communication  is  dominated  by  global  transposes  needed  for  the  3D 
FFTs.  This  involves  an  all-to-all  communication  pattern  that  minimizes  latency  and 
maximizes  message  size.  Such  communication  can  be  challenging  for  large  distributed 
memory  systems  and  must  be  coded  carefully.  It  also  requires  a  fast  inter-node  network. 

Extremely  high  parallel  efficiency  is  achieved  with  TRIPLE  by  minimizing  serial 
operations,  conducting  efficient  parallel  I/O  using  hierarchical  data  structuring  and  write- 
behind  data  buffering  (Werne,  et  al.  2000),  and  using  fast  communication  calls,  available 
either  via  SHMEM  or  MPI  library  routines.  Figure  1  shows  parallel  scaling  tests  using  a 
fixed  per-processor  problem  size  on  five  different  platforms.  The  variable  plotted,  C,  is 
proportional  to  the  wall-clock  time  per  compute  cycle.  Fluctuations  in  the  data  with 
NCPU  result  from  the  cache-reuse  characteristics  of  the  different  radix  FFT  routines  for 
different  problem  shapes.  In  the  left  panel,  C  is  multiplied  by  the  peak  FLOP  rate  S  of 
the  processor  used  so  that  the  communication  characteristics  of  the  different  platforms 
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can  be  isolated."  In  the  right  panel,  C  alone  is  plotted,  revealing  the  relative  wall-clock 
time  needed  to  complete  a  simulation  on  the  platforms  for  which  we  are  requesting 
Challenge  time  (P5+  and  XT3).  P4+  data  is  also  included  in  the  right  panel  so  that  the 
P4+  and  P5+  can  be  directly  compared.  Of  the  platforms  tested,  the  Cray  XT3  exhibits 
the  best  asymptotic  performance  with  100%  parallelization  efficiency  (i.e.,  C»S  is 
constant  within  the  fluctuations  in  the  data).  The  Cray  T3E,  IBM  P4+,  and  IBM  P5+ 
exhibit  similar  parallel  performance,  with  an  efficiency  of  99.976%  according  to 
Amdahl’s  Law  or  m=95%  when  evaluating  the  grind  time.2 3  The  P4+  exhibits  slightly 
better  parallel  performance  than  the  P5+  because  the  inter-node  network  for  the  two 
machines  is  basically  the  same,  while  the  P5+  has  twice  as  many  processors  per  node  and 
its  processor  is  1 1.8%  faster.4  The  Xeon  cluster  exhibits  the  worst  perfonnance  (99.89% 
Amdahl’s  Law  or  m=88%  for  the  grind  time). 

The  XT3’s  superior  performance  is  due  to  its  excellent  network  -  7.6  GB/s  inter¬ 
node  bandwidth  and  5. 4-6. 6  ps  latency.  The  T3E  and  P4+/P5+  have  bandwidths  of  420 
and  450  MB/s  and  latencies  of  5. 6-8. 2  and  2-3  ps,  respectively.  The  ratio  of  the  per- 
processor  bandwidth  to  the  per-processor  FLOP  rate  is  0.066,  0.35,  and  1.46  B/flop  for 
the  P4+,  T3E,  and  XT3,  revealing  why  the  XT3  is  advantageous  for  TRIPLE’s  all-to-all 
communication  pattern.  Nevertheless,  the  parallel  efficiency  for  the  T3E  and  the 
P4+/P5+  is  more  than  adequate  for  our  needs.  Once  the  faster  processor  speed  S  and 
heterogeneous  network  of  the  P5+  is  factored  in,  overall  the  P5+  is  only  10%  slower  than 
the  XT3  for  certain  “sweet  spot”  values  of  the  grid  size  and  NCPU  (e.g.,  NCPU=900, 
1200  and  1600;  see  right  panel).  We  take  advantage  of  this  behavior  when  conducting 
simulations  on  the  P5+  by  carefully  choosing  value  of  NCPU.  In  contrast  to  these 
systems,  the  Xeon  cluster  perfonnance  was  determined  to  be  unacceptable  for  TRIPLE. 
The  inter-node  bandwidth  per  processor  on  this  system  was  only  250  MB/s,  and  the 
latency  was  19  ps  when  both  processors  per  node  are  used.  As  a  result,  the 
computational  cost  grew  too  rapidly  with  NCPU  to  be  a  useful  platfonn  for  TRIPLE. 

3.1.2.  TRIPLE  Job-Preparation  and  Automation  Scripts 

Significant  logistical  problems  arise  from  many  of  the  supporting  operations 
associated  with  large-scale  computing.  We  have  addressed  these  challenges  by 
developing  a  modular  system  of  code  and  Perl  scripts  to  simplify,  streamline,  and 
automate  large-scale  computing  tasks.  Strategies  include:  1)  automating  job  preparation 
and  submission;  2)  defining  an  irreducible  set  of  job-related  parameters  for  the  job-input 
file  so  that  redundancy  is  removed  during  batch-job  specification  (and  user  error  thereby 
minimized);  3)  automatically  migrating  large  3D  data  volumes  to  archival  storage  in  real 
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"  If  no  time  were  spent  in  communication,  OS  would  equal  the  operation  count  per  mode 
per  time  step  per  processor,  and  all  data  in  the  left  panel  would  be  identical  (and  flat). 
Separate  curves  result  from  differences  in  the  way  communication  is  managed  by  the 
different  platforms. 

3  Grind  time  gn  is  defined  by  gn=gi/nm,  where  gi  is  the  time  required  for  a  calculation  on 
one  processor,  and  gn  is  the  time  used  on  n  processors.  Parallel  efficiency  is  given  by  m. 

4  On  average,  TRIPLE  is  3.6%  faster  on  the  P5+  than  on  the  P4+.  This  is  apparent  in  the 
right  panel  of  Figure  1  where  C  is  not  normalized  by  the  processor  speed  S. 
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time  so  that  impact  on  local  supercomputer  disk  space  is  minimized;  and  4)  developing  a 
translation  layer  for  predictable  tasks  so  that  code  can  be  more  easily  maintained  at 
multiple  centers  and  users  need  not  memorize  multiple  syntaxes  to  solve  the  same 
problem;  these  tasks  include  fde  transfer  operations,  queue-system  control  directives  and 
parameters,  and  job-submission  syntax.  This  last  set  of  tools  (item  4)  forms  the  core  of 
the  Uniform  Command  Line  Interface  (UCLI);  see  www.pstoolkit.org,  home  of  the 
Practical  Supercomputer  Toolkit.  With  this  suite  of  software  tools,  users  on  this  project 
were  able  to  avoid  many  repetitive  and  mundane  tasks,  and  instead  were  able  to 
concentrate  on  higher-level  scientific  tasks. 

3.1.3.  Technical  and  Computational  Challenges 

The  challenges  for  the  simulations  proposed  arise  from  the  anisotropic  and  episodic 
nature  of  stratospheric  turbulence  and  the  enormous  range  of  length  scales  it  possesses; 
see  Figure  2.  The  wind-shear  and  wave-breaking  processes  1)  develop  from  flow 
instabilities,  2)  progress  through  a  high-Re  laminar  state  before  3)  becoming  fully 
turbulent,  4)  subsiding  as  turbulence  decays,  and  5)  the  flow  restratifies  (Werne  et  al. 
2005).  Though  a  properly  specified  SGS  model  can  be  useful  in  mitigating  the  enormous 
cost  of  stage  3  in  this  scenario,  it  cannot  address  the  similarly  daunting  resolution 
requirements  that  accompany  stage  2,  and  though  it  might  be  able  to  handle  stage  4,  care 
must  be  taken  at  late  time  as  the  spectral  gap  closes  between  the  energy-containing  and 
dissipation  ranges  in  the  spectrum  and  the  SGS-model  assumptions  are  violated.  When 
the  SGS  model  loses  validity,  DNS  is  the  only  available  option,  and  it  demands  large 
computational  resources;  e.g.,  computer  memory  and  cost  scale  as  Re94  and  Re3, 
respectively.  This  is  prohibitive  when  one  considers  stratospheric  Reynolds  numbers  of 
10  -10  .  In  contrast,  current  DNS  of  stratified  turbulence  reach  values  of  only  Rc~  1 0  . 
Though  still  expensive,  thankfully  these  Reynolds  numbers  are  sufficient  (barely)  to 
access  the  asymptotic  range  of  turbulence  at  high  Re,  and  direct  comparisons  with 
atmospheric  observations  demonstrate  their  relevance  (Werne  and  Fritts  2000).  Until  we 
discover  how  to  exploit  ways  to  overcome  the  inherent  difficulties  presented  by  stage  2 
and,  to  a  lesser  degree,  stage  4,  expensive  DNS  calculations  will  continue  to  be  necessary 
when  studying  episodic  atmospheric  turbulent  phenomena. 


3.1.4.  DNS/LES  Data  Analysis 

The  DNS  and  LES  results  were  analyzed  in  number  of  ways.  First,  flow  visualization 
was  conducted  to  determine  the  basic  flow  morphology  and  dynamics.  Next,  an 
exhaustive  statistical  analysis  package  was  developed  to  quantify  turbulence  budgets  and 
other  quantities  of  interest.  These  included  spanwise  and  full-horizontal  averaged  terms 
in  the  evolution  of  all  of  the  Reynolds-stress  components  UiUj,  the  turbulent  heat  flux  Ui0, 
the  dissipation  tensors  Vui*Vuj  and  V0*V0,  as  well  as  other  quantities  (e.g.,  2nd-order 
structure  constants  for  all  flow  fields,  anisotropy  invariants  II  and  III,  and  more.  Here  Ui 
is  the  ith  component  of  the  fluctuating  velocity  field  and  0  is  the  fluctuating  potential 
temperature.  Results  were  compared  with  published  atmospheric  measurements  and  with 
data  from  field  campaigns  designed  to  look  at  stably  stratified  atmospheric  dynamics. 
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DNS  and  LES  results  were  compared  to  evaluate  the  skill  of  the  SGS  model  in 
reproducing  the  DNS  solutions. 

3.2.  Analysis  of  Observational  Data 

Observational  data  were  analyzed  to  better  understand  atmospheric  dynamics  that 
give  rise  to  optical  turbulence  in  the  upper  troposphere  and  lower  stratosphere  and  to 
validate  the  DNS  results.  High-resolution  rawinsonde  data  from  the  CASES-99  and 
VTMX  field  campaigns  were  used  to  evaluate  the  properties  of  turbulent  layers  observed 
in  the  atmosphere,  while  aircraft  data  collected  using  three  NOAA  BAT  probes  aboard 
the  Grob  G520T  Egrett  were  used  to  compare  in  detail  with  the  simulation  results.  The 
Egrett  data  was  also  used  to  deduce  shear-turbulence  characteristics  in  the  edges  of  jet 
streams  over  England  and  Australia. 

3.2.1.  Rawinsonde  Data 

The  rawinsonde  data  analyzed  were  collected  at  4  sites  (Smileyberg,  Leon, 
Beaumont,  and  Eldorado)  in  Kansas  during  October  1999  and  at  3  sites  (Wheeler  Fann, 
Jordan  Narrows,  and  the  Salt  Lake  City  Airport)  in  Utah’s  Great  Salt  Lake  Basin  during 
October  2000.  In  total,  data  from  440  sonde  flights  were  analyzed.  Measurement 
profdes  collected  included  atmospheric  pressure,  temperature,  relative  humidity,  GPS 
location,  and  wind  speed  and  direction  at  typically  a  few  meters  vertical  resolution.  The 
maximum  altitude  attained  varied  from  flight  to  flight,  with  the  highest  flights  reaching 
up  to  28  km. 

In  order  to  determine  observed  turbulent-layer  properties  from  the  rawinsonde  data, 
the  following  procedure  was  employed.  First,  quality  control  of  the  data  files  was  done, 
and  bad  data  were  identified.  Second,  where  possible,  linear  interpolation  was  used  to 
estimate  bad-data  values  for  altitude  and  pressure.  Third,  cubic-spline  interpolation  was 
used  to  interpolate  the  raw  data  to  a  unifonn  grid  with  a  vertical  resolution  set  to  the 
minimum  vertical  spacing  reported  in  the  raw  data  file;  the  resolution  was  reduced  if 
more  than  40,000  vertical  levels  resulted  from  this  procedure.  The  resulting  interpolated 
data  possessed  resolutions  of  order  O(lm).  Care  was  taken  to  address  the  finite  bit- 
resolution  of  the  measurements  by  identifying  instrument  threshold  values  for  over¬ 
sampled  data.  Without  this  step,  many  false  turbulent  layers  would  be  reported.  Vertical 
derivatives  were  then  taken  using  3 -point  Lagrangian  differentiation. 

From  the  processed  rawinsonde  data,  we  identified  atmospheric  shear  layers,  and  we 
computed  parameters  of  interest,  such  as  the  layer’s  altitude  Z,  depth  L,  velocity  jump 
AU,  potential  temperature  drop  AT,  turning-shear  ratio  Vo/Uo  (defined  below),  etc.,  and 
derived  quantities,  such  as  the  layer  Reynolds  number  ReL=AU  L/v  and  bulk  Richardson 
number  RiL=gaAT  L/(AU)  .  Here,  g,  v  and  a  are  the  acceleration  due  to  gravity,  the 
kinematic  viscosity,  and  the  thennal  expansion  coefficient  in  the  atmosphere.  Layers 
were  identified  by  first  locating  local  maxima  in  ozU*ozU,  where  U  is  the  horizontal  wind 
vector  and  dz  is  the  vertical  derivative.  Adjacent  local  minima  (above  and  below  the 
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computed  mid-layer  maxima)  or  zero-crossings  of  the  shear  in  the  direction  of  maximum 
mid-layer  shear  were  then  used  to  define  the  layer  edges. 

We  also  evaluated  the  degree  to  which  one  might  deduce  the  actual  atmospheric 
Richardson-number  profde  given  a  simulated  profde  predicted  by  a  model.  To  do  this  we 
needed  to  smooth  the  rawinsonde  data  to  successively  coarser  resolutions  before 
analyzing  the  results.  This  was  carried  out  on  the  interpolated  data  by  smoothing  to  1 1-, 
33-,  100-,  and  300-m  vertical  resolutions  using  box-car  averaging. 

The  Interactive  Data  Language  (IDL)  software  package  from  Research  Systems,  Inc., 
was  used  for  the  rawinsonde  data  analysis. 

3.2.2.  Aircraft  Data 

Aircraft  data  were  collected  by  Dr.  George  Hacker  and  colleagues  using  Airborne 
Research  Australia’s  high-altitude  Grob  G520T  Egrett.  Measurements  were  taken  at 
altitudes  up  to  15  km  at  a  mean  airspeed  of  100  m/s  and  with  a  maximum  flight  duration 
of  up  to  8  hours.  Three  NOAA  BAT  probes  were  located  under  the  wings  and  high  on 
the  tail.  They  reported  temperature  and  all  three  wind  components  at  50  Hz.  The  data 
were  analyzed  by  Drs.  Don  Wroblewski  (Boston  University)  and  Owen  Cote  (AFRL). 

Regions  of  interest  in  the  aircraft  data  were  initially  identified  by  noting  significant 
vertical  velocity  fluctuations  above  the  background  fluctuations  of  about  0.3  m/s.  The 
data  were  low-pass  filtered,  and  second-order  structure  constants  for  all  measured 
quantities  were  evaluated.  Also,  properties  of  cliff-ramp  structures  observed  in  the  data 
were  quantified,  and  results  from  the  measurements  were  compared  with  the  DNS  results 
described  above. 

By  examining  the  data  in  this  way,  only  a  portion  of  the  data  has  been  processed  and 
analyzed  to  date. 

3.3.  Probabilistic  Modeling 

An  unanticipated  development  at  the  outset  of  this  work  was  the  ability  to  combine 
our  DNS/LES  results  with  observational  data  to  produce  hybrid  deterministic/ 
probabilistic  Cn  and  SGS  models,  which  can  forecast  both  statistical  expectation  values 
and  model  uncertainty.  This  is  possible  because  of  progress  in  recent  years  in  applying 
Bayesian  Hierarchical  Modeling  (BHM)  techniques  to  geophysical  applications  where 
the  character  of  BHM  is  particularly  well  suited  to  addressing  measurement  and  model 
uncertainty,  while  quantitatively  incorporating  significant  scientific  understanding.  A 
review  and  examples  relevant  to  geophysics  are  given  in  Berliner  (2003).  Applications 
discussed  there  include  climate-change  and  space-time  modeling  of  near-surface  ocean 
winds;  also  see  Berliner  et  al.  2003. 

For  our  purposes  of  improving  the  ABL  ADA,  BHM  provides  opportunities  for 
hybrid  deterministic/probabilistic  modeling  and  forecasting  of  complex  phenomena  that 
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would  otherwise  be  extremely  difficult  or  impossible.  In  essence,  with  BHM, 
complicated  joint  probability  distributions  are  modeled  as  products  of  simpler, 
component  models.  The  strategy  separates  the  role  of  measurement  uncertainty  from 
process  variability,  and  it  offers  direct  opportunities  for  building  physics  directly  into  the 
component  models.  As  such,  BHM  is  a  natural  framework  for  organizing  the  various 
data  sources  needed  for  an  improved  ABL  ADA,  including,  e.g.,  rawinsonde, 
thermosonde,  and  aircraft  measurements,  radar,  satellite,  and  possibly  other  data, 
compiled  or  parameterized  results  from  high-resolution  DNS  and  LES  of  turbulence,  as 
well  as  theoretical  descriptions.  The  technique  is  equally  valid  for  1)  post-processing 
mesoscale-model  output  to  predict  Cn2  along  operational  paths  and  2)  as  a  basis  for  a 
probabilistic  SGS  scheme  inside  a  running  mesoscale  model. 

3.3.1.  Bayesian  Hierarchical  Modeling 

The  probabilistic  modeling  approach  we  employ  is  based  on  Bayes  Theorem 
(Bernardo  and  Smith  1994).  To  describe  the  method,  we  begin  with  the  basic  formula 

[A  |  F]  oc  J  [F  |  A,Y]  [A  |  Y]  [Y]  dY,  (2) 

where  [Y]  refers  to  the  probability  of  the  occurrence  of  event  Y,  [A|Y]  is  the  conditional 
probability  of  A  given  Y,  and  [F|A,Y]  is  the  joint  conditional  probability  of  F  given  A 
and  Y. 

The  variable  F  represents  the  relevant  subset  of  the  resolved  mesoscale-model 
variables,  such  as  the  forecast  temperature,  pressure,  wind  speed,  etc.;  i.e.,  F  is  our 
forecast  before  improvement  by  BHM.  The  variable  A  represents  the  set  of  SGS 
variables  being  predicted  by  the  model.  In  other  words,  A  represents  an  improvement  to 
the  forecast  F  via  BHM.  Finally,  the  variable  Y  represents  the  set  of  parameters  that 
characterize  the  underlying  (and  unresolved)  processes  that  give  rise  to  A;  e.g.,  if  the 
process  being  modeled  results  from  observed  atmospheric  turbulent  mixing  layers,  then  Y 
can  be  described  by  the  layer  depth  F,  velocity  jump  AU,  altitude  Z,  etc.  Alternately,  if 
the  process  involves  the  overturning  and  breaking  of  atmospheric  gravity  waves,  then  Y 
would  be  characterized  by  the  wave  amplitude,  frequency,  and  wavelength,  as  well  as  the 
background  stratification. 

The  conditional  probability  distribution  [A|Y]  appearing  in  (1)  is  called  the  physics 
prior  distribution  of  A  given  the  unknown  variables  Y.  The  conditional  joint  probability 
[F|A,Y]  in  (1)  is  the  so-called  likelihood  distribution.  Finally,  the  output  of  the  model 
[A|F]  is  called  the  posterior  distribution,  which  is  our  probabilistic  prediction  for  A. 

Applications  of  (2)  relevant  to  this  project  are  discussed  in  §4.3. 

4.  RESULTS  AND  DISCUSSION 

In  this  section,  we  summarize  results  for  our  DNS/FES  solutions,  atmospheric  data 
analysis,  and  the  new  BHM  approach  to  probabilistic  prediction. 
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4.1.  DNS/LES  Solutions 


The  numerical  solutions  analyzed  for  this  project  were  computed  over  the  past  ten 
years,  beginning  in  1997.  Early  validation  efforts  of  our  wind-shear  simulations  were 
conducted  with  smaller  numerical  domains  than  those  computed  most  recently. 
Nevertheless,  comparisons  with  published  and  more  recently  collected  atmospheric  data 
have  demonstrated  that  even  the  early  numerical  solutions  provide  a  quantitatively 
accurate  description  of  stably  stratified  atmospheric  dynamics.  In  §4.1.1.,  we  discuss 
results  from  these  early  validation  exercises,  and  we  explain  why  the  DNS  solutions  at 
lower  Reynolds  numbers  than  those  attained  in  the  upper  troposphere  and  lower 
stratosphere  are  relevant  and  useful  for  ADA  development.  After  addressing  solution 
validation,  we  next  describe  the  behavior  of  the  computed  wind-shear  solutions  in  more 
detail  in  §4.1.2.  In  §4.1.3.  and  §4.1.4.,  we  discuss  the  implications  of  these  results  for 
Cn2  and  SGS  turbulence  modeling. 

After  describing  our  wind-shear  solutions,  we  then  describe  our  gravity-wave 
breaking  simulations  in  §4.1.5.  (linear  stability  problem)  and  §4.1.6.  (nonlinear 
dynamics).  In  §4.1.7.,  we  discuss  the  implications  for  the  wave-breaking  solutions  on  Cn2 
and  wave-saturation  modeling. 

4.1.1.  Numerical-Solution  Validation 

Early  wind-shear  solutions  were  computed  in  domains  of  size  (A,,  A/3,  2X)  in  the 
(streamwise,  spanwise,  vertical)  directions,  where  A,  is  the  wavelength  of  the  most 
unstable  asymptotic  linear  eigen  mode  for  the  Kelvin-Helmholtz  instability.  Second- 
order  structure  functions  were  computed  for  the  temperature  T  and  all  velocity 
components  (U,V,W),  and  comparisons  were  made  with  available  atmospheric  data. 
Here  (U,V,W)  are  the  velocity  components  in  the  (x,y,z)  directions,  where  x  (y)  is 
aligned  in  the  streamwise  (spanwise)  direction,  and  z  points  in  the  vertical.  The  second- 
order  structure  function  A r\|/  for  the  variable  \| /  is  defined  by  Ar\j r  =  ((\|/(x+r)-\|/(x))~). 
Here,  r  is  the  spatial  separation  of  two  points  in  the  3D  field,  and  ( )  indicates  averaging 
over  the  spatial  domain  of  x.  Note,  that  for  very  small  spatial  separations  (i.e.,  r— >0), 
Ar\|/2  =  ((V\|/  •  t )2)  r2,  where  r  is  a  unit  vector  in  the  r  direction.  For  larger  spatial 
separations,  the  dependence  of  Arv|T  on  r  is  determined  by  the  spatial  statistics  of  \|/.  For 
example,  if  the  fluctuations  in  \|/  are  generated  by  homogeneous-isotropic  turbulence, 
then  Ar\|/  ~  r"  in  the  inertial  subrange  for  both  T  and  (U,V,W)  (Monin  and  Yaglom 
1975).  In  contrast,  for  strongly  stratified  turbulence,  Bolgiano  (1959)  predicts  ArT”  ~  r 
and  ArU2  ~  r6/5. 

We  evaluated  Ar\|/2  from  the  simulation  results  and  fit  to  the  form  Ar\|/2  =  Cv2  r“  for 
separations  r  in  the  r=x  and  r=y  directions  for  \\i  =  T,  U,  V,  and  W.  Here,  Cv  is  the 
second-order  structure  constant  for  \|/,  and  a  is  a  power-law  exponent.  We  compared  our 
fit  parameters  with  the  atmospheric  measurements  by  Kaimal  et  al.  (1972),  Gurvich  and 
Zubkovskii  (1966),  Paquin  and  Pond  (1971),  Wyngaard  and  Cote'  (1971),  Wyngaard  and 
Pao  (1971),  and  Gibson  (1962,  1963).  Since  none  of  these  atmospheric  measurements 
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evaluated  a  directly,  but  rather  assumed  a=2/3,  they  must  be  considered  as  calibrations 
for  the  coefficient  Cv  .  Nevertheless,  this  is  a  useful  comparison,  as  it  allows  a  measure 
of  how  well  the  simulation  results  capture  the  magnitude  of  the  inertial  sub-range  of 
atmospheric  turbulence. 

In  Werne  and  Fritts  (2000),  we  reported  that:  Ct2=3.3s'1/3x  for  both  r=x  and  r=y; 
Cu2=2.1s2/3  for  r=x  and  Cu2=(4/3)2.1s2/3  for  r=y;  Cv2=1.6s2/3  for  r=y  and  Cv2=(4/3)1.6s2/3 
for  r=x;  and  Cw  =(4/3)0. 9s  for  both  r=x  and  r=y,  during  times  when  the  dynamics  are 
nearly  statistically  stationary.  These  results  are  consistent  with  the  atmospheric  results 
for  Ct2,  which  are  Ci2=Ces'1/3x,  with  Ce=3.3±0.3  (Kaimal  et  al.  1972),  Ce«3.5  (Gurvich 
and  Zubkovskii  1966),  and  Ce»3.3  (Paquin  and  Pond  1971).  Similarly,  atmospheric 
measurements  of  Cu‘=Cs  also  agree  with  our  results,  reporting  C=2.0±0.1  (Kaimal  et 
al.  1972),  C=2.1±0.2  (Wyngaard  and  Cote'  1971),  C=2. 1+0.1  (Wyngaard  and  Pao  1971), 
and  CA2.1  (Gibson  1962,  1963).  Finally,  though  some  evidence  for  C=1.6  (for  V)  and 
C=0.9  (for  W)  was  previously  reported  in  the  literature  (Monin  and  Yaglom  1975),  these 
results  remained  somewhat  inconclusive  until  recent  aircraft  data  (Wroblewski  et  al. 
2003)  also  obtained  these  same  results.  It  is  important  to  note  that  this  observational 
verification  of  the  simulation  results  for  Cv~  and  Cw“  represents  much  more  than  a  simple 
calibration  of  the  numerical  results.  Indeed,  the  fact  that  C  for  U,  V,  and  W  are  not  all 
equal  (not  even  close)  clearly  demonstrates  that  atmospheric  turbulence  under  stably 
stratified  conditions  is  far  from  homogeneous  and  isotropic. 

Another  valuable  aspect  of  the  aircraft  measurements  reported  by  Wroblewski  et  al. 
(2003)  is  that  they  directly  evaluated  a,  and,  like  the  earlier  simulation  results  reported  by 
Werne  and  Fritts  (2000),  exponents  ranging  between  a=2/3  and  a=2/5  for  T  and  U  were 
also  observed.  It  should  be  noted  that  a=6/5  was  never  seen  in  either  the  simulations  or 
the  aircraft  data;  hence,  Bolgiano’s  1959  theory  clearly  does  not  apply. 

From  the  computed  2nd-order  structure  constants  and  scaling  exponents  a,  we  see  that 
the  DNS  solutions  adequately  describe  2nd-order  statistics  in  the  inertial  sub-range  of  the 
atmospheric  turbulence  spectrum  under  stable  conditions.  In  addition,  because  the  DNS 
does  not  model  small-scale  quantities  and  instead  must  resolve  all  relevant  physical 
scales  of  motion  (right  down  to  the  dissipation  scale),  we  are  also  able  to  evaluate  the 
spatial  inner  scale  i0  at  which  the  observed  a=2/3  or  a=2/5  laws  transition  to  the  a=2  law 
for  small  separation  separations.  We  did  this  for  all  fields  analyzed,  and  we  found  the 
following  for  T  and  U:  £0(T,r=x  or  r=y)  «  7.4  £K;  £0(U,r=x)  «  V  2  8.0  £K  and  £c(U,r=y)  « 

8.0  £k-  Here  /k=(v3/s)14  is  the  Kolmogorov  length  scale.  We  note  that  the  temperature 
inner  scale  f0(T)  was  measured  previously  in  the  atmosphere,  and  it  agrees  perfectly  with 
the  DNS  result  (Hill  and  Clifford  1978).  Also,  the  ratio  of  the  inner  scales  deduced  from 
the  longitudinal  £0(U,r=x)  and  transverse  £0(U,r=y)  2nd-order  structure  functions,  i.e., 

f 0(U ,r=x)/ f 0(U ,r=y)  *  V  2,  agrees  with  nearly  isotropic  motions  at  small  scales  because 
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<3yU0yU)  «  2(<9XU5XU),  and  hence  the  smallest  spatial  scale  in  x  for  U  is  larger  than  the 
smallest  spatial  scale  in  y  by  V  2. 

A  legitimate  and  potentially  perplexing  question  regarding  these  results  is  “How  can 
DNS  at  Re=3xl04  simultaneously  duplicate  2nd-order  structure  constants  and  the  inner 
scales  of  motion  (i.e.,  quantities  that  span  the  broad  range  of  length  scales  from  the  tiny 
dissipation  scale  up  to  the  largest  scales  of  turbulent  motion  in  the  inertial  sub-range) 
when  Re  can  be  as  large  as  Re=10  -10  in  the  atmosphere?”  The  answer  is  that  the 
simulations  have  achieved  sufficient  resolution  to  clearly  separate  the  smallest  and  largest 
length  scales  in  the  flow.  Once  such  scale  separation  is  achieved,  we  then  expect  the 
scaling  behavior  for  asymptotically  large  Re  to  become  apparent.  Hence,  flows  with  Re 
spanning  a  large  range  of  values,  but  which  are  all  in  the  same  asymptotic  scaling  regime, 
can  be  meaningfully  compared  once  the  power-law  exponents  (e.g.,  a)  and  the  amplitude 
coefficients  (e.g.,  C¥  )  are  known. 

4.1.2.  Wind-Shear  Smulations 

The  wind-shear  simulations  were  conducted  using  a  constant  background  temperature 
gradient  T=(3z  and  a  hyperbolic  tangent  initial  streamwise  velocity  profile 
U=Uotanh(z/h).  The  flow  was  initiated  with  the  most  unstable  asymptotic  linear  eigen 
mode  with  an  amplitude  equal  to  7%  of  the  maximum  initial  mean  vorticity  plus  a 
Kolmogorov  noise  field  with  a  vorticity  amplitude  of  1.4%.  Cases  have  been  run  with  a 
range  of  Re  and  Ri  values  and  in  domains  of  varying  sizes.  Most  recently,  we  have 
completed  runs  in  domains  of  size  (4k,  2k,  2k),  where  k  is  the  most  unstable  wavelength, 
for  the  following  Ri  and  Re  values:  (Ri,  Re)=(0.05,  2500),  (0.15,  2900),  and  (0.20,  4000). 
Here,  Re  was  chosen  so  that  all  solutions  required  nearly  the  same  number  of  spectral 
modes,  which  peaked  at  (nx,  ny,  nz)  =  (3600,  1800,  1800).  These  large-domain  solutions 
were  computed  in  order  to  accumulate  statistics  for  the  BHM  effort  described  in  §X.  We 
have  also  completed  an  exploratory  case  with  (Ri,  Re)=(0.10,  2500)  in  a  domain  of  (k,  k, 
2k)  so  that  we  may  determine  the  value  of  Re  to  use  in  the  larger  domain  (it  is  Re=2800), 
and  so  that  we  may  get  a  glimpse  of  the  solution  behavior  for  Ri=0. 10. 

Flow  Morphology  and  the  Transition  to  Turbulence 

The  overall  evolution  for  the  four  cases  discussed  above  is  summarized  in  Figures  3- 
10.  Figures  3-6  show  the  temperature  on  both  a  streamwise-vertical  cut  and  on  the  mid¬ 
plane  at  seven  distinct  times  during  the  flow  evolution.  Figures  7-10  show  the  vorticity 
magnitude  at  the  same  times  and  locations.  From  both  fields  (temperature  and  vorticity), 
we  see  the  clear  influence  of  stratification  as  Ri  is  varied.  Most  obvious  is  the  reduced 
layer  depth  for  higher  Ri.  This  is  because  the  stronger  stratification  at  higher  Ri  impedes 
vertical  motion,  suppressing  layer  vertical  growth.  In  addition  to  this  continuous  trend 
with  Ri,  we  also  see  that  a  fundamental  change  occurs  between  Ri=0.10  and  Ri=0.15. 
For  the  lower  values  of  Ri,  we  see  laminar  billow  cores  persist  later  in  the  flow  evolution, 
while  turbulence  erodes  the  braid  regions  between  billows  first.  Only  later  do  the  cores 
become  fully  turbulent.  In  contrast,  for  the  higher  values  of  Ri,  the  braids  persist  intact 
longer,  while  turbulence  develops  in  the  cores  first.  This  difference  has  potentially 
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important  consequences  for  both  mechanical  and  optical  turbulence  generation  and 
evolution. 

To  understand  how  this  shift  in  turbulence  initiation  between  the  billow  and  braid 
regions  can  have  important  consequences,  consider  that  persistent  billow  cores  act  as 
strong  organizing  agents  for  the  flow.  At  low  Ri,  the  billows  are  large  and  roughly 
circular  in  cross  section  and,  as  a  result,  they  rotate  more  like  solid  bodies  than  their 
higher  Ri  counterparts.  The  low-Ri  billows  therefore  are  stabilized  longer  by  a  flywheel 
effect,  and,  thus,  are  able  to  extract  more  kinetic  energy  longer  from  the  background 
shear  flow,  leading  to  more  intense  turbulence  and  vigorous  mixing,  which  impacts  the 
structure  and  evolution  of  the  layer’s  Cn2  profile. 

The  stabilization  of  the  billow  cores  by  rotation  also  leads  to  a  contorted  and 
somewhat  counter-intuitive  transition  to  turbulence  that  gives  rise  to  two  distinct  episodes 
of  intense  turbulence  eruption,  while  the  higher  Ri  cases  display  a  simpler  transition 
characterized  by  only  a  single  episode  of  intense  turbulence  growth  and  decay.  To 
explain  these  distinct  pathways  to  turbulence,  we  begin  with  the  simpler  case  of  high  Ri. 
In  this  case,  as  soon  as  billows  fonn  and  undergo  their  first  rotation,  the  initial 
background  stable  stratification  in  their  cores  becomes  convectively  unstable  as  the 
billows  rotate  through  180°.  Since  the  strong  stratification  at  high  Ri  suppresses  vertical 
motion  and  overturning,  these  high-Ri  billows  are  compressed  in  the  vertical  direction 
and,  hence,  they  are  not  stabilized  by  solid-body  rotation  like  their  lower-Ri  counterparts. 
As  a  result,  high-Ri  billows  develop  turbulence  in  their  cores  nearly  immediately  upon 
formation.  Meanwhile,  their  braid  regions  remain  intact.  This  is  intuitive,  as  the 
enhanced  vertical  compression  of  the  braid  region  (resulting  from  adjacent  billow 
formation)  makes  the  braids  the  most  (convectively)  stable  part  of  the  flow.  When  the 
turbulence  from  the  billow  cores  is  sheared  out  of  the  cores  and  migrates  into  the  braids, 
then  the  braids  themselves  become  unstable. 

In  contrast,  the  transition  to  turbulence  for  the  low-Ri  billows  is  more  complex. 
Though  the  billow  cores  are  the  most  convectively  unstable  regions  in  the  flow,  they 
remain  stable  due  to  their  solid  body  rotation.  This  allows  the  billows  to  wind  up  fluid  of 
alternating  density  and  vorticity  at  their  periphery  so  that  strongly  layered  density  and 
spanwise  vorticity  develops  while  the  billows  remain  laminar.  Eventually  the  local 
Rayleigh  number  in  the  billow  edges  becomes  sufficiently  high  for  convective  instability 
to  develop  in  the  edge  regions  of  the  billows.  This  instability  takes  the  form  of 
streamwise-aligned  counter-rotating  convection  roles  that  wrap  around  the  billows.  The 
billow  edges  therefore  become  turbulent  first,  triggering  instability  in  the  adjacent  braid 
regions,  despite  the  enhanced  stable  stratification  in  the  braids.  Billow-core  stability 
persists  because  of  the  stabilizing  influence  of  rotation,  but  eventually  the  billow  cores 
also  become  unstable,  and  when  they  do,  turbulent  kinetic  energy  spikes. 

These  contrasting  transitions  to  turbulence  are  evident  in  Figure  11,  which  shows  the 
time  evolution  of  the  volume-integrated  kinetic  (KE,  solid  curve)  and  an  estimate  of  the 
potential  energy  (PE,  dashed  curve)  per  unit  volume,  along  with  the  maximum  values  of 
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the  three  vorticity  components  co;  for  all  four  cases  discussed  above.  The  vertical  dotted 
lines  correspond  to  the  times  shown  in  Figures  3-10. 

In  Figure  11,  the  background  mean  velocity  and  temperature  have  been  removed, 
hence, 

KE  =  (LxLyh)  '  J  [(f/-f/)2  +  V2  +  W^x  (1) 

PE  =  (/..  A  /?)  Ri  J  (T- f  )Vx  (2) 

The  normalization  factor  LxLyh  includes  the  horizontal  extent  of  the  layer  LxLy  in  order  to 
correct  for  the  effect  domain  size  has  on  the  integrated  energy.  We  use  h  for  the  vertical 
length  instead  of  Lz  because  the  turbulent  motions  do  not  fill  the  full  depth  of  the 
simulated  domain.  The  vorticity  maxima  indicate  when  turbulent  fluctuations  and  mixing 
are  active. 

Panel  (a)  in  Figure  1 1  (for  Ri=0.05)  exhibits  the  two-peak  episodic  nature  of  KE  and 
PE  for  weak  stratification  described  above.  The  strong  second  peak  occurs  when 
turbulence  finally  erupts  in  the  billow  cores  for  this  case.  This  is  very  different  from  the 
Ri=0.10  case,  for  which  we  see  that  billow  and  braid  turbulence  occur  nearly 
simultaneously  (see  Figures  4  and  7)  ,  and  hence  only  one  large  peak  is  evident  in  both 
KE  and  PE.  For  Ri=0.15  and  0.20,  billow  turbulence  clearly  develops  first,  and  the 
braids  become  unstable  relatively  later  (see  Figures  5,  6,  8,  and  9).  In  these  cases,  there  is 
mild  evidence  for  the  distinct  timing  of  braid  and  billow  turbulence  in  PE,  but  the  effect 
is  not  nearly  as  pronounced  as  it  is  for  Ri=0.05.  The  reason  for  the  difference  is  that  the 
weakly  stratified  structure  of  the  billow  cores  is  conducive  to  the  production  of  copious 
amounts  of  turbulent  KE,  while  strongly  stratified  braids  are  not. 

Other  noteworthy  features  in  Figure  11  are  :  1)  the  coincidence  of  the  highest  small- 
scale  turbulence  intensity  (as  is  evident  in  the  largest  vorticity  fluctuations)  with  a  local 
suppression  of  PE  due  to  vigorous  turbulent  mixing;  2)  more  rapid  initial  KE  and  PE 
growth  for  lower  Ri;  3)  much  higher  KE  and  PE  maxima  for  low  Ri;  and,  curiously,  4)  a 
non-monotonic  dependence  of  the  peak  vorticity  on  Ri  (i.e.,  vorticity  fluctuations  peak  at 
values  of  approximately  70,  45,  60,  and  47  for  Ri  of  0.05,  0.10,  0.15,  and  0.20, 
respectively  -  see  Figure  11). 

Flow  Profiles  for  High  and  Low  Ri 

Additional  evidence  of  important  differences  in  the  flow  at  high  and  low  stratification 
can  be  seen  in  the  evolution  of  profiles  for  first-  and  second-order  properties.  These  are 
shown  in  Figures  12  and  13  for  Ri=0.05  and  Ri=0.20,  respectively.  The  left  panels  of  the 
figures  show  the  time  evolution  of  the  mean  temperature  deviation  T-z  and  the 
streamwise  velocity  as  functions  of  height,  z  is  the  initial  profile  for  T,  therefore  T-z  is 
the  deviation  of  T  from  its  initial  profile.  The  curves  for  Ri=0.05  (Figure  12)  and 
Ri=0.20  (Figure  13)  are  separated  by  20  and  10  time  units,  respectively.  The  buoyancy 
time  scale  for  these  two  cases  differs  by  a  factor  of  two,  explaining  the  difference  in  the 
time  units  we  chose.  The  middle  panels  in  the  figures  show  the  normalized  temperature 
and  velocity  variance,  i.e.,  the  variance  divided  by  its  maximum.  The  right  panel  shows 
the  maximum  used  to  normalize  the  variance  curves  in  the  middle  panels.  In  the  figures, 
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(qq)  =  <(U-  U(z))  +  V  +  W“)  is  the  sum  over  the  three  velocity  components.  Also, 
max(TT)  has  been  multiplied  by  Ri  in  the  right  panel  so  that  each  of  the  two  curves 
represents  a  contribution  to  the  total  energy,  which  is  approximated  by  E=(q2  +  Ri  T2)/2. 
As  is  evident  from  the  decay  with  time  shown  in  the  right  panel,  nonnalizing  the  variance 
curves  is  necessary  to  see  the  variation  with  height  of  the  profiles  at  late  times. 

The  figures  show  that  the  mean  velocity  and  temperature  are  modified  significantly 
during  the  course  of  the  flow  evolution.  The  deviation  in  the  mean  temperature  exhibits  a 
characteristic  S-shape,  with  the  Ri=0.05  (Ri=0.20)  layer  growing  from  2h  to  roughly  6h 
(3h)  in  depth;  i.e.,  the  higher  Ri  layer  experiences  significantly  less  growth.  The  mean 
velocity  in  both  cases  approaches  a  nearly  constant  gradient,  though  deviations  from  a 
strictly  linear  profile  at  late  times  are  easily  discernible,  especially  for  Ri=0.05.  This 
occurs  because  of  the  vigorous  mixing  at  midlayer  (see  the  lower  middle  panel  of  Figure 
12),  which  produces  positive  curvature  in  U(z)  at  z=0.  The  variance  in  temperature 
(upper  middle  panels)  exhibits  peaks  near  the  top  and  bottom  of  the  turbulent  layer.  This 
occurs  because  the  mean  temperature  gradient  is  expelled  from  the  middle  of  the  layer  by 
turbulent  mixing,  and  it  accumulates  at  the  edges  as  a  result.  Where  the  mean 
temperature  gradient  is  large,  so  too  is  the  potential  for  producing  temperature  variance. 
Note  the  inverse  relationship  between  the  temperature  and  velocity  variances  for  Ri=0.05: 
whereas  the  velocity  variance  peaks  near  midlayer  and  is  understandably  weaker  near  the 
layer  edges,  the  temperature  variance  has  the  opposite  behavior,  i.e.,  it  attains  its 
maximum  near  the  layer  edges,  and  is  minimal  near  the  middle  of  the  layer.  This  inverse 
relationship  is  notably  absent  in  the  higher  Ri  case  (Figure  13).  Instead,  for  this  case  (qq) 
and  (TT)  appear  to  have  nearly  identical  shapes.  Also  significant  is  the  much  smaller 
magnitude  for  (TT)  in  the  higher  Ri  case;  compare  the  right  panels  of  Figures  12  and  13. 
The  reasons  for  these  stark  differences  are  1)  the  relative  efficiency  of  mixing  for  the  two 
cases  (mixing  is  initially  much  more  efficient  for  lower  Ri)  and  2)  the  relative  depths  of 
the  two  layers  (lower  Ri  results  in  deeper  layers).  This  second  effect  insures  that  the 
higher  Ri  case  will  retain  strong  velocity  gradients,  especially  in  the  edge  regions  of  the 
flow  (see  lower  left  panel  of  Figure  13  at  the  middle  time  shown,  which  is  t=80).  As  a 
result,  the  local  stability  profile  (i.e.,  the  diagnosed  Richardson  number  versus  height)  for 
this  case  dips  back  below  the  stability  limit  of  Ri(z)=0.25  near  the  edges,  reinvigorating 
the  dynamics  in  the  edge  regions  and  giving  rise  to  renewed  turbulence  KE  production 
there.  This  is  remarkable  because  at  a  slightly  earlier  time  (t=70)  the  entire  stability 
profile  had  already  been  elevated  above  Ri=0.25,  indicating  that  the  time  of  active 
dynamics  had  already  passed. 

The  right  panels  of  Figures  12  and  13  clearly  demonstrate  that  the  lower  Ri  case 
(Ri=0.05)  experiences  extremely  rapid  instability  growth  and  higher  energies  compared 
to  the  high-Ri  case.5  We  also  note  that  stratification  has  a  much  more  evident  impact  on 
the  developing  dynamics  for  the  higher  Ri  case  because  PE  is  larger  than  KE  during  the 
initial  growth  phase  for  Ri=0.2;  this  is  not  the  case  for  Ri=0.05.  Finally,  the  decay  in 


5  In  these  panels,  the  time  scale  is  a  factor  of  two  larger  for  the  lower  Ri  because  the 
buoyancy  time  is  proportional  to  V  Ri. 
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energy  is  exponential  and  rapid  for  Ri=0.2,  while  it  is  much  more  gradual  for  Ri=0.05, 
and  it  appears  to  experience  a  shift  in  the  controlling  dynamics  at  t=250  (see  the  change 
in  the  decay  behavior  in  the  upper  right  panel  of  Figure  12). 

Universal  Aspects  of  Turbulent-Shear-Layer  Evolution 

Given  the  dramatic  differences  in  turbulent- shear-layer  evolution  for  strong  and  weak 
stratification  discussed  above,  the  prospect  of  developing  a  reliable  description  that  can 
be  used  as  the  bases  for  a  traditional  SGS  approach  is  daunting.  When  resolution  is 
coarse  (as  is  always  the  case  for  mesoscale  NWP  models),  the  SGS  modeling  challenge  is 
compounded  because,  in  this  case,  even  the  gross  behavior  of  the  shear  layer’s  large-scale 
dynamics  and  evolution  must  be  modeled;  i.e.,  the  modeling  task  is  not  simply  confined 
to  the  smallest  scales  of  motion  in  the  inertial  range  of  turbulence. 

Thankfully,  despite  the  significant  differences  that  accompany  different  Ri,  we  have 
discovered  through  our  DNS  work  that  certain  characteristics  of  shear-layer  evolution  are 
universal  or  nearly  universal.  To  demonstrate  this,  in  this  section,  we  describe  the  end- 
state  profiles  of  our  Ri=0.05  and  Ri=0.20  simulations  and  discuss  how  aspects  of  these 
profiles  can  be  correctly  anticipated.  To  begin,  Figure  14  shows  the  Ri(z)  profiles  for  the 
highest  and  lowest  Ri  simulations  (Ri=0.20  and  Ri=0.05).  The  dotted  (solid)  line  depicts 
the  initial  (final)  Ri  profile.  The  time  selected  for  the  final  profile  is  after  the  turbulent 
dynamics  have  subsided  and  the  flow  has  achieved  a  restratfied  end  state.  The  final  time 
for  Ri=0.05  is  a  factor  of  two  larger  than  that  for  Ri=0.20,  compensating  for  the  differing 
buoyancy  periods  for  the  two  cases;  i.e.,  when  time  is  scaled  by  the  buoyancy  frequency 
N,  the  end  time  Nt  is  identical  for  the  two  cases.  From  the  figure,  we  note  that  the  initial 
layer  is  only  1.29  times  deeper  for  the  Ri=0.05  case,  however,  by  the  final  time,  the  ratio 
of  the  layer  depths  is  significant;  i.e.,  the  Ri=0.05  layer  becomes  2.23  times  deeper  than 
the  Ri=0.20  layer  at  the  final  time.  The  fact  that  the  lower  Ri  case  results  in  a  deeper 
shear  layer  is  not  surprising,  given  that  the  stronger  background  stratification  for  the 
Ri=0.20  case  inhibits  deeper  growth. 

In  order  to  explore  the  evolution  of  Ri(z)  with  time  (not  just  the  initial  and  final 
states),  Figure  15a  shows  the  midlayer  (z=0)  Ri  values  for  these  two  cases.  Here,  we  see 
that  both  layers  seemingly  approach  similar  values  at  late  times.  It  is  interesting  that  this 
value  is  approximately  Ri=0.55,  i.e.,  well  above  the  neutral  stability  limit  of  0.25.  It  is 
also  somewhat  alarming  that  at  intermediate  times  the  evolutions  are  clearly  quite 
different,  and  in  particular,  that  the  Ri=0.20  case  significantly  overshoots  neutral  stability 
and  actually  rebounds  to  a  less  stable  end  state  than  it  attains  at  intermediate  times. 

In  an  attempt  to  explain  this  strange  behavior,  we  consulted  the  flow  visualizations 
(e.g.,  Figures  3-6)  and  detennined  that  the  Ri(z)  profile  is  inadequate  for  completely 
describing  the  layer  stratification.  The  difficulty  lies  in  Ri’s  definition,  which  includes 
only  vertical  gradients  of  temperature  and  velocity,  while  the  persistent  braids  for  the 
Ri=0.20  case  provide  significant  horizontal  strain-rate  gradients  that  can  feed  turbulence 
production  for  a  significant  amount  of  time.  Therefore,  we  considered  a  generalization  of 
Ri  by  defining  the  nondimensional  time-scale  ratio  N/S,  where  N  is  the  buoyancy 
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frequency,  and  S=(SySij)I/2  is  the  magnitude  of  the  large-scale  strain  rate  S1J=(5lUJ+6jUl)y. 
Here,  ( )y  refers  to  averaging  in  the  spanwise  direction  only.  We  note  that  because  of  the 
factor  of  V2  in  the  definition  of  S,j,  at  the  initial  and  final  times  when  the  mean  flow  is 
horizontally  homogeneous  and  gradients  in  z  are  the  only  gradients  present,  S  will  be 
only  half  the  size  of  3Z  U,  and  therefore  (N/S)  will  be  roughly  twice  as  large  as  Ri,  or 
N/S=V  2  Ri. 

Figure  15b  depicts  the  mid-layer  value  of  N/S  versus  Nt.  From  the  plot,  we 
immediately  notice  that  the  overshooting  evident  in  Figure  15a  has  been  removed, 
indicating  that  our  suspicion  that  neglecting  the  influence  of  horizontal  strain-rate 
gradients  was  indeed  to  blame  for  the  odd  behavior  of  Ri(z=0).  We  also  see  that  the  final 
restratified  value  for  N/S  appears  to  be  very  near  1.0.  This  is  extremely  encouraging,  as 
it  indicates  that  the  dynamical  evolution  of  the  layer  halts  when  the  time  scale  associated 
with  the  evolving  large-scale  strain  rate  (i.e.,  S'1)  drops  sufficiently  for  it  to  be  overtaken 
by  the  buoyancy  time  scale  (i.e.,  N'1).  Once  this  happens,  buoyancy’s  role  in  resisting 
vertical  motion  is  satisfied  and  a  balance  between  the  potential  destabilizing  influence  of 
large-scale  kinetic  energy  and  the  stabilizing  influence  of  stratification  is  achieved. 

Extending  this  analysis  further,  we  developed  a  very  simple  model  of  a  layer’s  final 
depth  based  solely  on  energetics  arguments.  If  one  assumes  horizontally  homogeneous 
initial  and  final  states,  it  is  trivial  to  show  that  the  final  layer  depth  L  should  scale  with 
N/S.  This  is  explicitly  demonstrated  with  Figure  16,  which  shows  L  versus  Ri  for  the 
DNS  solutions.  The  data  points  show  the  results  for  the  four  Ri  cases  considered,  while 
the  curve  is  a  fit  to  the  form  L=CiRi'1/2+  C2  with  Ci=1.5  and  C2—O.5. 

The  two  results  N/S«l  and  L~Ri'  will  be  useful  in  developing  the  BHM  model 
described  in  §3.3.1. 

4.1.3.  Gravity  Wave-Breaking  Simulations 

In  addition  to  the  wind-shear  simulations  discussed  above,  we  also  conducted  series 
of  wave-breaking  simulations  as  a  part  of  this  project.  The  results  from  the  wave¬ 
breaking  work  focused  on  instability  structures  and  the  transition  to  turbulence  resulting 
from  overturning  gravity  waves  and  the  role  such  breaking  waves  play  in  producing 
turbulence  layering  in  the  atmosphere  (Fritts  et  al.  2003,  2006,  2007). 

Figure  17  shows  a  sequence  of  images  of  vortex  tubes  that  develop  in  a  recently 
computed  wave-breaking  event.  The  left  column  of  images  show  the  side  view  of  the 
wave  with  the  numerical  domain  tilted  to  coincide  with  the  phase  of  the  wave,  while  the 
right  panels  show  the  end-view  of  the  edge  of  the  domain.  Though  the  wave  propagates 
relative  to  the  numerical  domain,  the  images  have  all  been  shifted  so  that  the  wave  phase 
is  shown  fixed,  with  the  primary  wave  sinusoidal  velocity  upward  and  to  the  left  in  the 
upper  half  domain  of  the  left  column  and  downward  and  to  the  right  in  the  lower  half 
(again,  describing  the  orientation  in  the  left  column).  The  amplitude  of  the  wave  A  is 
such  that  a  portion  of  its  phase  is  statically  unstable  by  10%;  i.e.,  A=l.l.  For  this 
solution,  instability  sets  in  first  just  above  the  region  of  lowest  static  stability  as  the  wave 
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propagates.  The  dominant  instability  structures  take  the  form  of  horizontally-aligned, 
streamwise-oriented  vortex  tubes,  similar  to  the  predicted  most-unstable  eigen  mode 
computed  for  this  particular  wave  (Lombard  and  Riley  1996).  Also,  significant 
interactions  with  secondary  waves  are  also  observed.  The  initial  vortex  tubes  interact  and 
trigger  a  cascade  to  smaller  scale  motions,  rapidly  erupting  into  very  intense  turbulence  at 
small  scales.  The  turbulence  moves  through  the  flow  as  the  wave  propagates. 

Many  similarities  exist  between  the  wind-shear  results  discussed  above  and  the  wave¬ 
breaking  solutions  described  here;  these  include:  1)  the  presence  of  streamwise  aligned 
vortex  tubes  which  derive  KE  from  buoyancy  and  shear  sources,  2)  transitions  to 
turbulence  via  vortex  interactions  common  to  both  flows,  3)  expansion,  mixing,  and 
influences  of  turbulence  far  beyond  the  initial  regions  of  instability,  4)  turbulence 
duration  of  several  buoyancy  periods  or  less,  5)  sharp  gradients  in  turbulence  quantities  in 
edge  regions,  and  6)  spatial  separation  of  the  peak  turbulent  kinetic  energy  and  thermal 
dissipation  (indicating  separation  of  strongly  mixed  regions  from  Cn2  concentrations). 
This  last  effect  implies  that  optical  turbulence  impacts  will  not  be  coincident  with  the 
most  intense  regions  of  mechanical  turbulence,  but  instead  will  adjoin  such  regions. 
Important  differences  between  the  wind-shear  and  wave-breaking  results  are  also  present, 
the  most  immediate  being  that  wind-shear  events  are  localized  in  altitude,  while  waves 
propagate.  As  a  result,  wind-shear  turbulence  permits  the  development  of  sharper  edge 
regions  and  larger  peaked  concentrations  of  Cn2;  while  low-frequency  gravity  waves  may 
also  have  this  characteristic,  higher  frequency  waves  propagating  at  shallower  angles 
move  continuously  through  altitude,  sweeping  turbulence  through  the  fluid  following  the 
phase  of  the  breaking  wave. 

The  most  important  result  from  this  work  is  the  very  similar  evolution  of  a  companion 
simulation  conducted  with  A=0.9  (see  Figure  18).  Though  clear  differences  between  the 
two  solutions  exist  (Fritts  et  al.  2007),  conventional  wave-saturation  theory  (Dunkerton 
1989)  predicts  that  this  second  simulation  should  not  exhibit  instability  at  all.  This  is 
because  wave-saturation  modeling  takes  the  view  that  waves  that  are  everywhere 
statically  stable  are  also  dynamically  stable.  The  work  of  Lombard  and  Riley  (1996)  has 
relatively  recently  shown  that  this  idea  is  errant;  however,  they  also  demonstrate  that  the 
linear  growth  rates  for  low  amplitude  waves  is  lower  than  high  amplitude  waves,  so  it  has 
been  unclear  until  now  (with  our  simulations)  just  how  misleading  current  linear 
saturating  modeling  is.  Most  importantly,  while  saturation  modeling  would  predict  that 
the  A=0.9  wave  would  remain  stable  while  the  A=l.l  wave  would  experience  an 
amplitude  reduction  to  A=1.0,  and  then  become  stable,  our  simulations  show  that  actually 
both  solution  experience  significant  amplitude  reductions  to  A=0.3  before  saturating. 
This  is  very  far  from  conventional  saturation  models,  and  indicates  how  important 
simulations  like  these  will  be  in  developing  more  reliable  and  accurate  models  in  the 
future.  For  more  details  on  the  wave-breaking  solutions,  see  Fritts  et  al.  (2007). 

4.1.4.  Implications  of  DNS  Simulations  for  Cn2  and  SGS  Modeling 

Figure  19  shows  a  sequence  of  thermal  and  viscous  dissipation  images  for  wind-shear 
simulations  with  Ri=0.05.  This  image  clearly  communicates  the  nature  of  mixing  in 


18 


stratified  flows.  Whereas  mixing  is  vigorous  where  the  viscous  dissipation  is  large 
(orange  regions),  the  thermal  dissipation  is  weak  in  these  regions.  This  is  because  strong 
gradients  in  velocity  have  homogenized  the  thermal  field  at  the  core  of  the  billow,  and 
thermal  gradients  therefore  are  relegated  to  the  edges  of  the  layer.  Since  thermal  and 
index-of-refraction  gradients  are  intimately  related,  this  result  reveals  that  optical 
turbulence  will  also  be  concentrated  in  the  regions  shown  in  blue  in  the  figure. 

This  is  problematic  for  optical  turbulence  forecasting  because  the  edge  regions  of  a 
stratified  shear  layer  are  particularly  challenging  to  model.  First,  these  regions  are 
situated  in  the  transition  zone  between  quiescent  laminar  flow  exterior  to  the  shear  layer 
and  turbulent  flow  inside  the  layer.  As  a  result,  the  typical  mathematical  assumptions 
made  in  turbulence  modeling  (namely  homogeneity,  isotropy,  and  stationarity)  are  all 
known  to  be  false  in  this  region.  Furthermore,  stratification  remains  challenging  to 
model,  even  at  small  scales,  and  the  edges  of  the  layer  harbor  the  highest  mean 
stratification  levels  at  intermediate  and  late  times.  This  means  a  sophisticated  approach 
to  modeling  is  required  in  this  region.  We  note  that  the  new  results  shown  in  Figure  13  at 
high  stratification  levels  (i.e.,  Ri=0.20)  indicate  that  strong  velocity  fluctuations  also 
reside  in  the  edge  regions  at  intermediate  times.  This  potentially  amplifies  the  challenge 
because  it  implies  that  accurate  solutions  require  realistic  descriptions  of  both  the 
temperature  and  velocity  fields,  and  one  cannot  make  a  weak-flow  assumption  for  the 
velocity,  which  might  have  been  possible  at  lower  Ri  (see  Figure  12). 

4.1.5.  LES  Solutions  and  Turbulence  Modeling. 

As  part  of  this  project,  we  also  conducted  LES  of  the  same  problems  for  which  we 
computed  DNS  solutions.  Our  initial  plan  was  to  refine  the  LES  SGS  model  sufficiently 
for  stratified  flow  to  use  the  LES  code  in  place  of  the  much  more  expensive  DNS  code  to 
compute  solutions  and  accumulate  the  statistics  needed  for  the  BHM  parameterization 
approach.  After  some  refinement  of  the  LES  code  and  detailed  comparisons  with  the 
DNS  solutions,  we  made  progress  on  describing  the  dynamics  during  the  periods  when 
the  turbulence  is  most  intense,  but  had  difficulty  getting  the  solutions  to  match  when  the 
LES  was  started  from  t=0  with  the  same  initial  conditions  as  the  DNS. 

Figure  20  demonstrates  this  by  presenting  results  for  three  sets  of  comparisons.  The 
top  row  of  panels  in  the  figure  shows  results  for  a  TKE  SGS  model  for  the  LES.  The  left 
panel  shows  the  resolved-scale  kinetic  energy  (KE)  for  the  DNS  (red)  and  LES  (blue) 
solutions.  The  right  panel  shows  the  eddy  viscosity  (red)  and  eddy  diffusivity  (blue)  for 
the  LES  (solid  line)  and  filtered  DNS  (data  symbols)  fields.  Though  the  LES  appears  to 
track  the  DNS  resolved-scale  KE  reasonably  well  on  average  (upper  left  panel),  it  is 
clearly  much  too  viscous  and  diffusive,  with  eddy  transport  coefficients  that  are  8.5  times 
larger  than  the  DNS  values  (upper  right  panel).  The  dynamic-Smagorinsky-model  results 
(middle  row  of  panels)  are  far  superior  at  estimating  the  eddy-transport  coefficients, 
being  only  10-15%  too  large  when  turbulence  is  most  intense  (from  t=120  to  180),  and 
only  a  factor  of  2  too  large  (on  average)  at  earlier  and  later  times.  Nevertheless,  the 
dynamic-Smagorinsky  model  clearly  under-estimates  the  resolved-scale  KE  when  3D 
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instability  develops  (from  t=50  to  t=80)  and  at  later  times  during  the  turbulent  decay 
phase  of  the  flow  (middle-left  panel). 

Some  causes  for  concern  for  the  LES  approach  include:  1)  the  initial  deviation  of  the 
solutions  just  before  t=50,  when  3D  instability  of  the  flow  first  appears,  2)  the  subsequent 
detachment  of  the  phase  of  KE  oscillations  for  the  DNS  and  LES  solutions  (middle-left 
and  top-left  panels  from  t=80  to  t=  1 50),  3)  the  over-estimate  of  SGS  turbulent  transport  at 
late  times,  and  3)  large  spikes  in  the  LES  SGS  transport  coefficients  that  do  not  exist  in 
the  DNS  solutions.  Though  some  of  these  effects  appear  to  be  systematic  and 
representative  of  the  SGS  models  used,  others  may  be  a  direct  result  of  the  inadequacy  of 
the  modest  domain  used  (A  x  a/3  x  2 A)  to  produce  sufficiently  stable  statistics  so  that 
meaningful  DNS/LES  comparisons  can  be  made.  For  this  reason,  and  to  ensure  that  we 
will  be  able  to  use  the  DNS  to  accumulate  statistics  for  the  BHM  approach,  we  conducted 
simulations  in  much  larger  domains  of  size  4A  x  2A  x  2 A,  i.e.,  a  full  24  times  larger. 
Results  for  the  DNS  and  the  dynamic-Smagorinsky  LES  are  shown  in  the  bottom  row  of 
Figure  20. 

From  Figure  20,  we  see  that  KE  oscillations  are  significantly  reduced  in  the  large 
domain  (bottom-left  panel),  indicating  the  significantly  improved  statistics  for  the 
increased  domain  size.  The  departures  of  the  LES  from  the  DNS  that  remain  can  be 
clearly  identified,  and  speculation  regarding  the  potential  role  of  inadequate  spatial 
statistics  is  no  longer  a  concern.  The  differences  between  LES  and  DNS  that  persist 
include  1)  a  systematic  under-estimate  of  resolved-scale  KE  from  just  after  the  onset  of 
3D  instability  (t=37)  to  the  eruption  of  turbulence  in  the  billow  cores  (t=86),  2)  a 
systematic  over-estimate  of  resolved-scale  KE  from  t=  110  onward  during  the  turbulence 
decay  phase  of  the  evolution,  and  3)  sizeable  spurious  spikes  in  the  LES  eddy-transport 
coefficients  during  the  development  of  3D  instability. 

The  overestimate  of  the  eddy  coefficients  when  3D  instability  develops  is 
understandable  and  is  actually  impossible  for  the  LES  to  avoid.  It  results  when  sharp 
laminar  gradients  form  that  simply  must  be  resolved  by  increased  grid  resolution.  Unable 
to  resolve  these  high-Re  laminar  features,  the  SGS  model  interprets  them  as  resulting 
from  turbulent  motions,  and  eddy  transport  (erroneously)  spikes  as  a  result.  Furthermore, 
with  eddy  transport  grossly  over  estimated,  KE  dissipation  is  elevated,  and  resolve-scale 
KE  is  errantly  depleted  as  a  result.  Therefore,  the  LES’  inability  to  properly  describe  the 
onset  and  development  of  the  3D  instability  structures  in  the  flow  produces  both  the  large 
spike  in  the  eddy-transport  coefficients  and  the  dip  in  resolved-scale  KE.  It  seems  clear 
that  no  improvement  to  the  LES  SGS  model  will  be  able  to  address  this  difficulty,  and 
DNS  will  be  required  to  represent  this  period  in  the  flow’s  evolution.  This  is  unfortunate, 
as  the  resolution  requirements  during  this  time  are  as  high  as  at  any  later  time,  even  when 
turbulence  is  at  its  most  intense. 

The  overestimate  of  resolved-scale  KE  during  the  turbulent  decay  phase  is  interesting 
because  it  is  not  evident  for  either  of  the  LES  SGS  models  in  the  smaller  domain.  This 
indicates  that  this  aspect  of  the  numerical  solutions  cannot  be  reliably  evaluated  and 
analyzed  in  the  smaller  domain.  Continued  analysis  of  the  new  solutions  in  the  larger 
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domain  is  required  to  understand  the  over  estimate  of  resolved-scale  KE  (and  PE  -  not 
shown)  and  the  simultaneous  over  estimate  of  the  eddy-transport  coefficients  at  late  time. 
Differences  between  these  solutions  for  Ri=0.05  and  the  other  Ri  values  must  also  be 
evaluated. 

Despite  these  deficiencies  with  the  LES  solutions,  there  are  reasons  for 
encouragement.  Most  importantly,  the  differences  in  second-order  quantities  (like 
resolved-scale  KE)  between  the  DNS  and  LES  solutions  are  modest.  Hence,  we  expect 
second-order  structure  functions  to  be  reasonably  well  described  by  the  LES.  Indeed, 
Figure  21  presents  comparisons  of  fit  parameters  for  the  second-order  structure  functions 
for  the  DNS  and  LES,  and  the  results  appear  to  agree  sufficiently  for  the  LES  results  to 
be  directly  useful.  This  is  an  exciting  result,  as  the  LES  employed  6  times  fewer  grid 
points  in  each  spatial  direction  compared  to  the  DNS,  translating  into  a  factor  of  1300 
less  computer  time  required  to  obtain  the  result.  Nevertheless,  optimism  should  be 
tempered  until  an  observed  dependence  on  the  filter  width  is  understood  and  can  be 
anticipated  ahead  of  time.  Until  then,  calibration  of  the  LES  (via  judicious  choice  of  the 
number  of  grid  points  used)  will  be  required. 


4.2.  Atmospheric  Data  Analysis 

The  BHM  approach  described  in  §3.3.  and  §4.3.  requires  statistical  information  that 
meaningfully  quantifies  1)  the  likelihood  of  turbulence-generating  processes  in  the 
atmosphere  and  2)  the  nature  of  the  turbulence  resulting  from  these  processes.  The 
results  of  the  data  analysis  presented  in  this  section  directly  address  the  first  point,  by 
measuring  probability  density  functions  (PDFs)  that  characterize  atmospheric  turbulence 
layers,  and  indirectly  addresses  the  second,  by  providing  comparisons  with  the  DNS/LES 
solutions  so  that  we  may  be  confident  the  simulated  dynamical  processes  are  directly 
relevant  to  atmospheric  dynamics. 

Two  data  sources  were  used  in  the  data  analysis  associated  with  this  project.  These 
include  1)  rawinsonde  balloon  data  and  2)  aircraft  measurements  of  atmospheric 
turbulence.  The  data  and  analysis  methods  are  described  in  §3.2.  Here,  we  present  the 
results  of  those  analyses. 

4.2.1.  Balloon-Data  Analysis 

Figure  22  shows  probability  density  functions  resulting  from  the  analysis  procedure 
described  in  §3.2.1.  The  upper  left  panel  shows  the  PDF  for  the  layer  depth  L,  which  was 
obtained  by  locating  the  local  maxima  and  minima  of  the  strain  rates  <3ZU  and  <3ZV 
associated  with  the  measured  wind  field  (see  §3.2.1.).  The  upper  middle  and  lower  left 
panels  show  similar  results  for  the  velocity  jump  AU  and  the  layer's  altitude  Z.  The 
remaining  three  panels  show  additional  information,  including  the  temperature  drop  AT 
(upper  right  panel),  the  layer's  bulk  Richardson  number  Rfi  (lower  right  panel),  and  the 
layer  Reynolds  number  ReL  (lower  middle  panel).  Each  plot  shows  three  curves,  where 
the  solid  curve  represents  results  obtained  when  using  all  of  the  data,  while  the  dashed 
and  dotted  curves  are  for  the  stratosphere  and  troposphere  only. 
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These  initial  results  are  extremely  encouraging  for  using  probabilistic  methods  to 
evaluate  atmospheric  quantities.  In  particular,  all  of  the  PDFs  are  well-defined,  and  in 
addition  those  for  L,  AU,  and  AT  have  nearly  perfect  universal  shapes.  To  see  this,  in 
Figure  23,  we  have  re-plotted  the  results  using  the  reduced  variable  X’=(X-(X))/XRMs, 
where  X  represents  any  of  the  variables  plotted,  and  (X)and  Xrms  are  its  mean  and  root- 
mean-squared  values,  respectively.  When  the  PDF  is  invariant  in  its  reduced  variable  X’, 
it  can  be  parameterized  solely  by  its  mean  and  RMS  values  (X)  and  XRMs,  significantly 
simplifying  the  modeling  effort.  Nevertheless,  the  mean  and  RMS  values  differ  notably 
between  the  troposphere  and  stratosphere,  and  they  exhibit  some  small  variation  with 
location  and  probably  season  as  well.  For  example,  stratospheric  layers  are  less  likely  to 
be  deep,  but  they  have  higher  velocity  and  temperature  drops  across  them.  Their 
Reynolds  numbers  are  less  than  for  the  troposphere,  consistent  with  the  troposphere's 
weaker  stratification  (and  therefore  its  more  vigorous  turbulence).  More  work  is  required 
to  determine  and  parameterize  the  seasonal  and  locational  dependencies.  However,  the 
remarkable  similarities  we  observe  when  data  from  individual  sites  (from  Utah  to 
Kansas!)  are  compared,  indicates  that  the  dynamics  aloft  are  fairly  robust  and  should  be 
successfully  modeled  by  this  kind  of  approach. 

4.2.2.  Aircraft-Data  Analysis 

Figure  24  shows  a  plot  of  cliff-ramp  structures  in  the  measured  potential  temperature 
reported  by  Wroblewski  et  al.  (2003)  for  a  flight  over  Wales  through  the  bottom  of  the  jet 
stream  at  an  altitude  of  1 1 .4km.  The  measured  atmospheric  data  are  shown  in  the  left 
column  of  panels  for  (from  top-to-bottom)  potential  temperature  0,  streamwise  velocity 
U,  spanwise  velocity  V,  and  vertical  velocity  W.  Distinctive  structures  like  these  have 
been  observed  in  1 1  data  sets  collected  over  Australia  and  England,  and  more  continue  to 
be  discovered  as  analysis  of  the  data  progresses. 

The  form  of  the  structures  shows  a  clear  sawtooth  shape  in  0  that  coincides  with  a 
subsequent  increase  in  U  and  suppression  of  fluctuations  in  W.  The  sawtooth  pattern  in  0 
is  immediately  identifiable  from  Figures  5  and  6  as  the  temperature  jump  in  the  braid 
region  between  billows.  As  the  braids  remain  intact  longer  for  the  higher  Ri  cases,  it  is 
more  natural  to  identify  this  sawtooth  feature  with  the  Ri=0.15  and  Ri=0.20  cases  we 
have  computed.  The  right  column  of  panels  in  Figure  24  shows  that  this  is  indeed  the 
case  for  a  path  through  the  DNS  data  at  t=77.8  h/Uo  for  Ri=0.20.  Similar  features  are 
evident  for  the  Ri=0.15  case  at  t»72  h/Uo,  but  such  events  become  more  rare  in  the  DNS 
as  Ri  is  decreased  further,  and  they  are  completely  absent  by  Ri=0.05.  Similarly,  the 
suppression  of  W  fluctuations  near  the  braid  region  is  a  characteristic  of  high-Ri  solutions 
where  instability  sets  in  first  in  the  billow  cores  (see  Figures  9  and  10).  In  contrast,  the 
lower  Ri  cases  possess  fluctuations  that  are  substantial  near  the  braids,  and  at  Ri=0.05  a 
large-scale  variation  in  W  is  observed.  This,  combined  with  the  absence  of  cliff-ramp  or 
sawtooth  structures,  rules  out  Ri=0.05  as  a  candidate  value  of  Ri  for  the  aircraft  data  in 
Figure  24. 
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We  can  refine  these  qualitative  observations  further  in  an  attempt  to  assign  a  specific 
value  of  Ri  to  the  observations.  We  do  this  by  quantitatively  comparing  more  features  in 
the  flight  data  to  match,  in  fingerprint  fashion,  similar  paths  through  the  DNS.  We  begin 
by  deducing  the  wavelength  X  evident  in  the  observations.  By  evaluating  data  from  the 
three  BAT  probes,  it  was  determined  that  the  aircraft  traveled  roughly  12°  off  of  the 
streamwise  direction  (Wroblewski  2007).  With  a  flight  speed  of  97.5  m/s,  the  time-trace 
data  is  converted  to  spatial  data,  and  X  is  detennined  to  be  A,=1.63  km.  Also,  from  the 
temperature  measurements  during  ascent,  the  potential  temperature  gradient  is  measured 
to  be  P=7.0  K/km.  Finally,  from  the  first  sawtooth  at  x=74.6  km  in  Figure  24,  we  obtain 
the  following  discernable  characteristics:  a  temperature  drop  of  AT=2.8  K,  a  velocity 
jump  of  AU=5.5  m/s,  and  a  drop  in  AT  by  0.75  at  x=79.6  km  (corresponding  to  the 
passage  of  At=51.3  seconds). 

Similarly,  from  the  DNS  results  in  Figure  24  for  Ri=0.20  (and  from  similar  results  for 
Ri=0.15)  [and  for  Ri=0. 10],  we  obtain  the  following:  X=9.5  h  (10.3  h)  [11.3  h];  AT=0.9 
ph  (1.65  ph)  [0.6  ph];  AU=0.64  U0  (0.7  U0)  [0.3  U0];  all  at  t=77.8  h/U0  (72.2  h/U0)  [81.2 
h/Uo];  and  a  reduction  in  AT  to  75%  of  its  initial  value  over  a  time  span  of  At=5.2  h/Uo 
(3.9  h/Uo)  [3.7  h/Uo].  While  these  quoted  quantities  were  straightforward  to  detennine 
from  the  virtual  flight  paths  through  the  DNS  for  Ri=0.20  and  Ri=0.15,  they  were  more 
challenging  to  obtain  from  the  Ri=0.10  case  (and  admittedly  required  some  imagination). 

From  these  data,  we  can  estimate  the  value  of  h  for  the  observed  shear  layer  by 
setting  the  DNS  and  aircraft  wavelengths  equal  to  each  other.  In  this  way  we  obtain 
h=172  m  (158  m)  [144  m]  from  the  Ri=0.20  (0.15)  [0.10]  DNS  results.  Similarly,  by 
equating  the  DNS  and  aircraft  temperature  drops  AT,  we  obtain  P=18.1K/km  (10.7K/km) 
[32.4K/km],  From  the  direct  measurement  of  P=7. OK/km,  we  see  that  Ri  must  be 
between  Ri=0.10  and  Ri=0.15.  Furthennore,  we  note  that  with  the  dramatic  increase  in 
the  prediction  of  p  for  the  Ri=0. 10  case,  we  conclude  that  Ri  cannot  be  as  low  as  Ri=0. 10 
for  the  aircraft  measurements.  By  simply  extrapolating  the  Ri=0.15  and  Ri=0.20  cases, 
we  obtain  Ri=0.13  for  the  aircraft  data.  Also,  if  we  perfonn  similar  analysis  attempting 
to  match  the  decay  of  AT  over  time,  we  estimate  Ri  must  be  greater  than  Ri=0.1 1. 

More  analysis  and  comparisons  between  the  DNS  results  and  the  aircraft 
measurements  are  required  to  refine  these  estimates  further.  This  is  part  of  our  ongoing 
work,  because  comparisons  like  these  will  permit  us  to  construct  the  census  infonnation 
needed  for  the  BHM  approach  described  in  the  next  section. 

4.2.3.  BHM  Approach  to  Probabilistic  Prediction 

The  basic  approach  of  a  BHM  was  described  briefly  in  §3.3.1  with  Equation  (2), 
which  contained  the  variables  A,  F,  and  Y.  In  the  approach,  A  represents  the  quantity  we 
are  modeling.  For  this  project,  it  is  a  measure  of  optical  turbulence  (e.g.,  A=Cn  )  or  the 
SGS  quantities  needed  to  improve  the  AFWA  mesoscale  model  (e.g.,  A=T;j  or  h;,  where 
Tq  is  the  SGS  Reynolds  stress  and  h;  is  the  SGS  heat  flux).  F  is  the  set  of  mesoscale- 
model  forecast  variables  (e.g.,  pressure,  temperature,  wind  fields,  etc.)  or  diagnostic 
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quantities  based  on  these  (e.g.,  the  model-predicted  Cn  profile,  or  (Cn")M,  the  model- 
predicted  Ri  or  RiM,  etc.).  Lastly,  Y  describes  the  unresolved  processes  that  contribute  to 
the  existence  of  A.  For  example,  for  atmospheric  mixing  layers,  Y  could  include  the 
layer  depth  L,  its  altitude  Z,  velocity  jump  and  temperature  drop  AU  and  AT,  initial  layer 
Richardson  number  Ri(0),  the  layer  Reynolds  number  ReL,  the  background  stratification 
N“,  and  the  layer’s  age  a.  For  gravity  waves,  Y  would  encompass  the  spectrum  of 
wavelengths  and  frequencies  present,  their  associated  amplitudes,  and  the  wind  and 
stability  profiles  through  which  the  waves  will  propagate. 


When  considering  these  lists  of  variables  proposed  to  define  Y,  immediately  we  note 
that  some  of  the  quantities  are  redundant.  For  example,  in  the  case  of  atmospheric 
mixing  layers,  if  we  know  the  location  Z  of  a  layer,  then  it’s  mean  viscosity  v(Z)  can  be 
estimated;  hence,  knowing  AU,  L,  and  Z,  we  can  compute  the  layer  Reynolds  number 
ReL  =  AU  L  /  v(Z),  demonstrating  that  ReL  is  not  a  linearly  independent  variable  required 
in  the  set  Y  if  Y  already  contains  AU,  L,  and  Z,  Furthermore,  if  we  know  ReL  and  Ri(<)), 
then  we  can  use  the  DNS  solutions  to  determine  the  non-dimensional  layer  depth  L/h  and 
the  non-dimensional  temperature  drop  AT/(oz  T  h)  as  a  function  of  time.  Using  these 
quantities  and  the  fact  that  (N/S)~  «  Rif  ~  Vi  (where  Rif  is  the  final  layer  Richardson 
number  -  see  §4.1.2  and  Figure  15b  to  understand  the  value  «  Vi),  then  we  can  estimate  h 
(from  L/h)  and  N2  from  the  approximation 
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which  follows  directly  from  Rif  «  Vi.  Noting  these  relationships,  we  can  remove  all 
redundancy  in  the  proposed  set  Y  and  use  only  the  five  variables  Y=(AU,  L,  Z,  Ri(0),  a). 
The  detailed  implementation  of  this  set  for  Y  is  accomplished  via  repeated  application  of 
Bayes  Theorem,  so  that  evaluation  of  [Y]  =  [L,  AU,  Ri(0),  a,  Z]  can  be  organized  by 
computing  the  individual  component  distributions  contained  in 

[Y]  =  [L  |  AU,Ri(0),a,Z]  [AU  |  Ri(0),a,Z]  [Ri(0)  |  a,Z]  [a|Z]  [Z]  (4) 
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The  final  step  in  developing  a  BHM  for  Cn  involves  detennining  which  set  of 
forecast  variables  F  in  Equation  (2)  are  sufficient  to  construct  the  algorithm.  Decisions  as 
to  which  variables  are  needed  will  be  based  on  the  results  of  experimentation  once  the 
model  is  implemented  and  tested.  The  set  F  will  evolve  as  we  learn  more  from  the  results 
of  these  tests. 


Even  though  we  do  not  know  which  variables  F  will  be  included  in  the  final 
implementation  of  (2),  we  can  still  describe  the  procedure  for  building  the  conditional 
PDFs  identified  in  (2),  and  we  can  describe  the  procedure  for  solving  the  model. 

First,  the  physics  prior  distribution  [Y]  has  already  been  detennined  in  the  limit  of 
Equation  (4)  when  all  component  variables  in  Y  are  independent,  i.e.,  Equation  (4) 
reduces  to  [Y]=[L][AU][Ri(0)][a][Z].  The  individual  component  distributions  for  this  case 
are  presented  in  Figures  22  and  23.  More  detailed  statistical  analysis  of  the  rawinsonde 
data  is  required  to  evaluate  Equation  (4)  more  completely.  Next,  the  physics  prior 
distributions  [Cn2| Y]  must  be  computed  from  the  DNS.  Preliminary  evaluations  have 
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already  been  computed  based  on  second-order  structure  functions  from  the  DNS  (Weme 
and  Fritts  2000).  As  was  stated  in  §4.1.1.,  these  DNS  second-order  structure  function 
results  have  been  validated  against  atmospheric  tower  data  and  aircraft  data,  and  they 
have  been  shown  to  be  in  perfect  agreement  with  the  atmospheric  measurements. 
Finally,  the  remaining  likelihood  distribution  [  F  |  Cn2,  Y  ]  on  the  right-hand  side  of 
Equation  (2)  represents  an  estimate  of  the  dependence  of  the  relevant  forecast  variables  F 
on  the  actual  Cn2  field  and  the  Y  associated  with  unresolved  atmospheric  mixing  layers. 
This  can  easily  be  determined  from  compiled  rawinsonde  data  by  1)  detennining  the 
properties  Y  from  the  data  (which  we  have  done  when  constructing  the  PDFs  presented  in 
Figures  22  and  23)  and  then  simply  filtering  the  raw  data  to  the  mesoscale  model’s  spatial 
resolution,  and  computing  F  from  the  result.  For  example,  if  a  method  for  estimating  the 
model-predicted  Cn“  profile,  i.e.,  (Cn  )m,  is  used  [e.g.,  Dewan  et  al.  (1993)  or  Jackson 
(2004),  models  currently  used  by  the  Air  Force],  then  F=(Cn2)M,  and  computation  of  [  F  I 
Cn  ,  Y  ]  is  nothing  more  than  a  filtering  procedure  applied  to  (Cn  )m- 

4.2.4.  Intrinsic  Limitations  of  Deterministic  NWP  Modeling 

Mesoscale  NWP  forecast  models  rely  on  SGS  descriptions  of  unresolved  processes. 
These  models  correctly  acknowledge  that  increased  atmospheric  stability  leads  to  reduced 
small-scale  vertical  transport  and  mixing.  As  such,  the  SGS  models  used  all  include  a 
dependence  on  Ri  for  the  model  mixing  length,  and  they  effectively  turn  off  SGS 
turbulent  transport  when  the  atmospheric  Richardson  number  is  above  that  needed  for 
instability  processes  to  become  active.  For  the  gradient  Richardson  number,  this  would 
result  for  flows  with  Ri>0.25.  However,  the  models  typically  use  somewhat  larger  values 
of  the  bulk  Richardson  number  to  control  the  SGS  model  dynamics. 

Because  of  the  importance  of  Ri  in  determining  the  SGS  dynamics,  a  BHM  based  on 
Equation  (2)  for  SGS  turbulence  will  contain  the  likelihood  distribution  [RiM  |  Ri],  and, 
therefore,  the  detenninistic  skill  of  the  method  will  hinge  on  the  degree  to  which  the 
model-predicted  Richardson  number  RiM  is  an  appropriate  proxy  for  Ri.  Furthermore, 
because  Ri  can  be  constructed  in  a  straightforward  manner  from  the  rawinsonde  data  (by 
smoothing  the  measured  wind  and  potential  temperature  fields  to  the  model  resolution, 
then  evaluating  the  derivatives  needed  to  compute  RiM),  we  can  determine  a  priori  when 
mesoscale  NWP  model  output  is  most  likely  to  produce  detenninistic  predictability  for 
turbulence  forecasts  and  when  it  is  not. 

Figure  25  shows  a  sample  of  the  rawinsonde  data  collected  during  the  DOE  Vertical 
Transport  and  Mixing  (VTMX)  field  program  conducted  in  Utah’s  Great  Salt  Lake  Basin. 
The  left  panel  shows  the  potential  temperature  9(z)  and  filtered  Ri(z)  profiles  (which  we 
denote  Rig(z))  for  a  single  balloon  launch.  Ri§(z)  is  plotted  on  the  right  for  three  different 
filter  widths  8.  The  shaded  regions  in  the  left  panel  indicate  altitudes  where  Rig(z)<0.25, 
with  larger  8  indicated  by  darker  shading.  The  shaded  regions  identify  turbulent  layers, 
with  darker  layers  likely  more  intense.  Because  Rig(z)  systematically  increases  with  8  as 
a  direct  result  of  reduced  gradient  resolution  by  the  coarser  grid,  the  higher-S  profiles 
predict  many  fewer  active  mixing  layers  than  are  actually  present.  This  is  significant 
because  the  coarsest  filter  width  used  8=300m,  typical  of  mesoscale-model  vertical 
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resolutions,  explicitly  demonstrates  that  NWP  models  based  on  current  common  practices 
alone  cannot  possibly  predict  but  a  small  fraction  of  the  dynamical  layers  that  naturally 
occur  in  stratified  environments;  hence,  it  is  the  case  that  conventional  NWP  model 
output  must  be  integrated  with  other  data  sources  to  exhibit  predictive  power  when 
attempting  to  forecast  atmospheric  turbulence. 

The  right  panels  in  Figure  25  show  normalized  PDFs  for  Ri,  i.e.,  [Ri]  for  the 
troposphere  and  stratosphere  combined.  The  top  panel  shows  the  results  from  a  single 
balloon  launch,  with  four  sets  of  data  indicated  for  each  of  four  filter  widths.  The  data  in 
the  bottom  panel  results  after  combining  data  for  82  separate  balloon  flights;  note  the 
much  smoother  profdes  and  clear  systematic  dependence  on  8  for  the  larger  statistical 
sample.  As  with  the  other  PDFs  examined,  the  normalized  Ri  PDFs  exhibit  a  universal 
shape  and  can  therefore  be  characterized  solely  by  their  mean  and  RMS  values. 

Figure  26  examines  the  conditional  probability  density  [  RiM  |  Ri  ],  estimated  from  the 
rawinsonde  measurements,  by  presenting  the  mean  and  RMS  of  the  ensemble  PDF’s 
presented  in  the  lower- right  panel  of  Figure  25.  Figure  26  shows  results  for  the  same 
filter  widths  8  used  in  constructing  Figure  25.  It  is  apparent  that  as  8  is  increased,  the 
relationship  between  the  mean  and  RMS  predicted  values  for  Rig  and  the  actual  value  Ri 
becomes  weaker,  i.e.,  the  curves  become  flatter.  When  the  curves  are  completely  flat, 
fidelity  from  the  procedure  to  predict  mean  and  RMS  values  is  lost,  and  the  estimated 
[RiM  |  Ri]  distribution  becomes  identical  to  that  predicted  by  a  regional  climatology 
model  for  Ri  based  on  archived  rawinsonde  data.  We  note  that  this  appears  to  occur  at  all 
of  the  rawinsonde  launch  sites  at  VTMX  and  CASES-99  at  a  resolution  of  about  300m. 
Hence,  in  order  to  obtain  high-fidelity  Ri  predictions  from  model  data,  NWP  models 
must  achieve  greater  resolution  than  they  do  currently,  and  the  evaluated  RiM  measures 
currently  used  in  stability-corrected  mixing-length  applications  undoubtedly  relate  poorly 
to  the  actual  values  that  should  be  used.  Clearly  a  BHM  approach  would  significantly 
improve  the  situation,  and  given  the  apparent  universal  form  of  the  ensemble  PDFs  in  the 
lower-right  panel  of  Figure  25,  Ri  climatologies  appear  to  be  feasible.  This  is  excellent 
news  for  BHM  methods  applied  to  regions  where  sparse  or  no  current  data  is  available, 
but  for  which  compilations  of  representative  PDFs  are  possible,  either  through  available 
archived  results  or  via  results  for  a  dynamically  similar  surrogate  region. 
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Figure  1.  Parallel  performance  on  five  supercomputer  platforms.  C  is  proportional  to  the  wall 
time  per  compute  cycle;  S  is  the  peak  processor  speed.  Left  panel  shows  C*S ,  which  isolates  the 
inter-processor  network  performance.  Asymptotic  fits  to  the  XT3,  P4+  and  Xeon  data  are  shown 
(dotted  lines).  SHMEM  or  MPI  communication  libraries  are  used,  depending  on  which  is  fastest. 
Right  panel  shows  C  alone,  which  indicates  the  relative  wall-clock  time  to  complete  an  operation 
on  three  of  the  platforms.  For  certain  NCPU  ( e.g .,  900,  1200  and  1600),  the  P5+  and  XT3  times 
are  within  -10%  of  each  other. 


Figure  2.  Vorticity  in  stratified  wind  shear  for  Ri=0.15  showing  large  range  of  length  scales 
present.  The  ratio  of  sizes  for  the  largest  and  smallest  eddies  is  roughly  150.  Given  that  5  to  8 
grid  points  are  required  to  resolve  the  smallest  eddies,  it  is  clear  that  very  large  numerical 
meshes  are  required  to  simulate  such  flows  with  DNS.  Only  two  of  the  four  large-scale  billows 
simulated  are  shown  so  that  small-scale  features  can  be  viewed.  Layer  Reynolds  numbers  of  Re  > 
30,000  are  required  to  achieve  sufficient  scale  separation  to  be  atmospherically  relevant,  while 
four  billows  are  needed  to  obtain  adequate  statistics  to  populate  the  probability  density  functions 
needed  for  the  BHM  approach  employed.. 


27 


a 


28 


Figure  3.  Temperature  in  stratified  wind  shear  for  Ri=0.05  and  Re=2500  at  seven  distinct  times 
during  the  flow  evolution.  For  each  time,  two  panels  are  shown:  1)  a  side  view  of  the  four 
billows  and  2)  a  top  view  of  40%  of  the  mid-plane.  The  times  chosen  correspond  to  when  (a) 
billows  reach  maximum  amplitude,  t=37;  (b)  secondary  instability  triggers  the  transition  to 
turbulence,  t=54;  (c)  KE  and  PE  dip  to  local  minima  before  secondary  peaks  in  KE  and  PE 
occur,  t=68;  (d)  KE  and  PE  exhibit  secondary  maxima,  t=85;  (e)  turbulence  intensity  and 
vorticity  magnitude  attain  maximum  values,  coincident  with  a  local  minimum  in  PE,  t=lll,  (f) 
turbulence  first  becomes  horizon  tally  homogeneous,  t-130;  and  (g)  turbulence  decay  exhibits  a 
change  in  character,  t=251.  Panels  a-d  appear  on  the  preceding  page. 
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Figure  4.  (Continued  next  page.) 
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Figure  4.  Temperature  in  stratified  wind  shear  for  Ri=0.10  and  Re=2500.  Seven  times  are 
presented,  as  in  Figure  3,  but  because  the  evolution  depends  strongly  on  Ri,  some  of  the 
times  correspond  to  different  stages  in  the  flow  evolution  and  morphology.  Times  chosen 
correspond  to  when  (a)  billows  reach  maximum  amplitude,  t=44;  (b)  secondary  instability 
appears  in  the  billow  edge  region,  t=54;  (c)  turbulence  intensifies  in  billow  edges  and 
reaches  braids,  t=68;  (d)  turbulence  peaks  in  the  billow  cores,  t=76;  (e)  vorticity 
magnitude  attains  maximum,  coincident  with  a  local  minimum  in  PE,  t=89,  (f)  turbulence 
first  becomes  horizontally  homogeneous,  t=97;  and  (g)  turbulence  decay  exhibits  a  change 
in  character,  t=165.  Panels  a-d  appear  on  the  preceding  page.  Note:  this  solution  is  only  of 
a  single  billow  which  has  been  replicated  to  be  consistent  with  Figures  3,  5,  and  6. 
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Figure  5.  Temperature  in  stratified  wind  shear  for  Ri=0.15  and  Re=2900.  Seven  times  are 
presented.  Times  chosen  correspond  to  when  (a)  billows  reach  maximum  laminar  amplitude, 
t=44;  (b)  billows  reach  maximum  amplitude,  with  instability  developing  in  the  billow  edge 
region,  also  coincident  with  maximum  PE,  t=57;  (c)  turbulence  invades  billow  cores,  t=67; 
(d)  turbulence  reaches  braids,  t-77;  (e)  turbulence  begins  exponential  decay  phase,  t=86; 

(f)  turbulence  first  becomes  horizontally  homogeneous,  t=91;  and  (g)  turbulence  decay 
exhibits  a  change  in  character,  t=116.  Panels  a-d  appear  on  the  preceding  page. 
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Figure  6.  (Continued  next  page.) 
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Figure  6.  Temperature  in  stratified  wind  shear  for  Ri=0.20  and  Re=4000.  Seven  times  are 
presented.  Times  chosen  correspond  to  when  (a)  billows  reach  maximum  laminar 
amplitude,  t=37;  (b)  billows  reach  maximum  amplitude,  with  turbulence  erupting  in  the 
billow  cores,  t=54;  (c)  PE  peaks,  t=66;  (d)  turbulence  reaches  braids,  t=82;  (e)  turbulence 
begins  exponential  decay,  t=90;  (f)  turbulence  first  becomes  horizontally  homogeneous, 
t=105;  and  (g)  turbulence  late  in  the  exponential  decay  phase,  t=130.  Panels  a-d  appear  on 
the  preceding  page. 
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Figure  7.  Vorticity  magnitude  in  stratified  wind  shear  for  Ri=0. 05  and  Re=2500.  Panels  a-d 
appear  on  the  preceding  page.  See  Figure  3 ’s  caption  for  times  associated  with  each  panel. 
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Figure  8.  (Continued  next  page.) 
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Figure  8.  Vorticity  magnitude  in  stratified  wind  shear  for  Ri=0.10  and  Re=2500.  Panels  a-d 
appear  on  the  preceding  page.  See  Figure  4  ’s  caption  for  times  associated  with  each  panel. 
Note:  this  solution  is  only  of  a  single  billow  that  has  been  replicated  to  be  consistent  with 
Figures  7,  9,  and  10. 
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Figure  9.  Vorticity  magnitude  in  stratified  wind  shear  for  Ri=0.15  and  Re=2900.  Panels  a-d 
appear  on  the  preceding  page.  See  Figure  5 's  caption  for  times  associated  with  each  panel. 
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Figure  10.  Vorticity  magnitude  in  stratified  wind  shear  for  Ri=0.20  and  Re=4000.  Panels 
a-d  appear  on  the  preceding  page.  See  Figure  6 ’s  caption  for  times  associated  with  each 
panel. 
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Figure  1 1  Time  evolution  of fluctuation  kinetic  energy  KE  (solid  line)  and  an  estimate  for  the 
potential  energy’  PE  (dashed  line)  (top  panels)  and  the  maximum  of  each  vorticity  component 
(bottom  panels)  for  (a)  Ri=0.05,  (b)  Ri=0.10,  (c)  Ri=0.15,  and  (d)  Ri=0.20.  Vertical  dotted 
lines  correspond  to  the  times  shown  in  Figures  3-10. 
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Figure  12.  Ri=0.05  -  Horizontally  averaged  profiles  for  the  mean  temperature  (upper  left)  and  mean 
streamwise  velocity  (lower  left)  and  for  the  normalized  temperature  variance  (upper  middle)  and 
normalized  velocity  variance  (lower  middle).  The  variance  normalization  factors  used  in  the  middle 
panels  are  shown  in  the  rightmost  panel.  The  solid  (dashed)  curve  depicts  the  velocity  variance 
(t q~ )=(u  +v~+w  )  (temperature  variance  (r ). 
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Figure  13.  Ri=0.20  -  Horizontally  averaged  profiles  for  the  mean  temperature  (upper  left)  and  mean 
streamwise  velocity  (lower  left)  and  for  the  normalized  temperature  variance  (upper  middle)  and 
normalized  velocity  variance  (lower  middle).  The  variance  normalization  factors  used  in  the  middle 
panels  are  shown  in  the  rightmost  panel.  The  solid  (dashed)  curve  depicts  the  velocity  variance 
(q2 )=  (u2 +v2 +w2 )  (temperature  variance  (T2 '). 
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Figure  14.  Ri(z)  profiles  for  wind-shear  simulations  with  Ri=0.20  (left)  and  Ri=0.05  (right). 
Initial  (dotted)  and  final  (solid)  profiles  are  shown.  The  time  of  the  final  profile  is  selected  at 
late  time  during  final  restratification. 
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Figure  15.  Mid-layer  Richardson  number  Ri(O)  (left panel)  and  N/S  (right panel)  versus 
Nt  for  Ri=0. 05  (dashed)  and  Ri=0.20  (solid).  Ri(O)  results  from  averaging  Ri(z)  from  z=- 
0.5  to  z=0.5.  See  text  for  defin  ition  of  S. 
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Figure  16.  Final  layer  depth  L  versus  initial  layer  Richardson  number  Ri  for  four 

wind-shear  simulations.  Data  points  show  the  simulation  results,  while  the  solid 

1  /2 

curve  is  a  fit  of  the  form  L=1.5  Ri~ -  0.05. 
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Figure  17.  3D  evolution  of  vortex  tubes  (X2  visualization)  for  an  amplitude  A= 1.1 
gravity  wave  viewed  from  the  side  and  from  above.  Red  and  yellow  regions 
denote  the  most  intense  rotation.  Wave  propagation  is  to  the  left  in  all  panels. 
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Figure  18.  3D  evolution  of  vortex  tubes  (A 2  visualization)  for  an  amplitude  A= 0.9 
gravity  wave  viewed  from  the  side  and  from  above.  Red  and  yellow  regions 
denote  the  most  intense  rotation.  Wave  propagation  is  to  the  left  in  all  panels. 
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Figure  19.  3D  volume  rendering  of  the  thermal  (blue)  and  viscous  (orange)  dissipation  for  a 
wind-shear  simulation  with  Ri=0. 05.  The  top  row  shows  the  view  from  above  at  three 
distinct  times.  The  bottom  row  shows  the  view  from  the  side. 
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Figure  20.  DNS/LES  comparisons  of  wind-shear  simulations  with  Ri=0.05  and  Re=2500. 
Left  panels  show  KE  for  LES  (blue)  and  DNS  (red).  Right  panels  show  eddy  viscosity  (red) 
and  eddy  diffusivity  (blue)  for  the  LES  (solid  line)  and  DNS  (data  symbols)  solutions.  Eddy 
quantities  are  evaluated from  the  DNS  by  first  filtering  to  the  LES  resolution,  and  then 
computing  the  SGS  quantities  directly  from  the  sub-filter  fields.  The  top  row  shows  results 
for  an  LES  with  a  TKE  SGS  model,  while  the  middle  and  bottom  panels  show  results  for  a 
dynamic-Smagorinsky  model.  The  top  and  middle  rows  show  results  for  DNS  and  LES  on 
small  domains  of  size  (X,  k/3,  2k)  using  up  to  (800,  270,  1600)  spectral  modes  for  the  DNS 
and  (120,  40,  240)  grid  points  for  the  LES,  while  the  bottom  row  shows  results  on  a  large 
domain  of  size  (4  k,  2  A,  2  k)  with  up  to  (3000,  1500,  1500)  modes  for  the  DNS  and  (480,  240, 
240)  modes  with  the  LES. 
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Figure  2 1  Structure-function-fit  parameters  versus  height  for  LES  (blue)  and  DNS  (red) 
solutions.  Structure  functions  are  computed  separately  for  streamwise  (x)  and  spanwise  (y) 
spatial  separations.  The  two  solutions  are  statistically  indistinguishable ,  demonstrating  that 
the  LES  is  able  to  reproduce  the  DNS  results  while  using  1300  times  less  computer  time.  See 
Werne  &  Fritts  2000 for  details  on  the  computation  of  the  2nd-order  structure-function-fit 
parameters. 
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Figure  22.  Probability  density  functions  for  atmospheric  shear  layers  derived from  440 
rawinsonde  flights  over  Kansas  and  Utah.  Top  row  fro  left  to  right:  L,  layer  depth ;  AU, 
velocity  jump  across  layer;  AT,  temperature  drop  across  layer.  Bottom  row  from  left  to 
right:  Z,  layer  height;  Re^AUL/v,  layer  Reynolds  number;  RiL=ga  AT/(AU)2,  layer  bulk 
Richardson  number.  All  panels  contain  three  curves.  The  dashed  (dotted)  curve  is  for  raw 
data  located  above  (below)  the  tropopause.  The  solid  curve  results  from  all  the  data. 
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Figure  23 .  Same  as  Figure  22,  but  all  PDFs  are  with  respect  to  the  reduced  variables 
X’=(X-(X))/Xrms,  where  X  represents  one  of  the  variables  plotted  in  Figure  22,  and  (X)and 
Xms  are  the  average  and  RMS  values  of  X  respectively.  If  the  PDF  versus  the  reduced 
variable  X’  is  identical  in  the  troposphere  and  stratosphere,  then  it  possesses  a  universal 
form  and  can  be  easily  represented  by  (X)  and  Xms- 
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Figure  23 .  Same  as  Figure  22,  but  all  PDFs  are  with  respect  to  the  reduced  variables 
X’=(X-(X))/Xrms,  where  X  represents  one  of  the  variables  plotted  in  Figure  22,  and  (X)and 
Xms  are  the  average  and  RMS  values  of  X  respectively.  If  the  PDF  versus  the  reduced 
variable  X’  is  identical  in  the  troposphere  and  stratosphere,  then  it  possesses  a  universal 
form  and  can  be  easily  represented  by  (X)  and  Xms- 
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Figure  24.  Cliff-ramp  structures  in  aircraft  measurements  (left)  and  DNS  wind-shear  solutions 
(right).  The  left  panel  shows  four  time  traces  (from  top  to  bottom)  for  aircraft-measured  potential 
temperature  6,  streamwise  velocity  U,  spanwise  velocity  V,  and  vertical  velocity  W  along  a 
streamwise  flight  track  over  Wales  at  and  altitude  of  11.4  km.  The  data  were  collected  using 
Airborne  Research  Australia,  ’s  Grob  G520T  Egrett  outfitted  with  three  NOAA  BAT  probes 
positioned  under  each  wing  and  high  on  the  tail  (see  Wroblewski  et  al,  2003).  Similar  structures 
are  observed  along  paths  through  the  DNS  at  Ri=0. 15  and  Ri=0.20.  The  right  panels  show  such 
structures  at  t=7  7  h/Uo-  See  text  for  more  discussion  of  the  data  and  comparisons  between  the 
DNS  and  aircraft  data 
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Figure  25.  Left  panel:  VTMX  Wheeler  Farm  rawinsonde  data,  17  Oct  2000,  0:48  IJTC.  Left-most  sub¬ 
panel  shows  raw  potential  temperature  0(z).  Rifz)  is  plotted  on  right  for  three  different  data  filters: 

8 -33m,  101m,  and  300m.  Rifz)  is  computed  after  applying  filter  to  G(z)  and  winds.  Regions  with 
Ri<0.25  are  shaded  gray,  with  the  shading  level  determined  by  the  value  of  the  largest  filter  scale  for 
which  Ri  <  0.25  is  satisfied.  Upper-right  panel:  Probability  Density  Functions  (PDFs)  of  Rifz)  data 
from  the  left  panel.  Solid  symbols:  1 1  m-smoothed  data.  Open  circles:  33m  smoothing.  Triangles: 
101m  smoothing.  Squares:  300m  smoothing.  Lower-right  panel:  Ensemble  PDFs  for  all  Wheeler 
Farm  data  collected  at  VTMX.  Note  the  reduced  error  bars  and  better-defined  functional  dependence 
on  5 for  the  ensemble  data. 
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Figure  26.  Mean  and  RMS  of  the  ensemble  PDF's  shown  in  Figure  25. 
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